Harnessing elasticity to generate self-oscillation via an electrohydrodynamic instability

TL;DR

This paper demonstrates how elastic structures in dielectric fluids can induce self-oscillations through electrohydrodynamic instabilities, offering new insights for designing adaptive soft machines.

Contribution

It introduces a novel elasto-electro-hydrodynamic model showing how elasticity transforms bifurcations into self-oscillations, supported by numerical and theoretical analysis.

Findings

Elastic structures can induce self-oscillations via EEH instability.

The model predicts critical conditions matching numerical results.

Elasticity transforms pitchfork bifurcation into Hopf bifurcation.

Abstract

Under a steady DC electric field of sufficient strength, a weakly conducting dielectric sphere in a dielectric solvent with higher conductivity can undergo spontaneous spinning (Quincke rotation) through a pitchfork bifurcation. We design an object composed of a dielectric sphere and an elastic filament. By solving an elasto-electro-hydrodynamic (EEH) problem numerically, we uncover an EEH instability exhibiting diverse dynamic responses. Varying the bending stiffness of the filament, the composite object displays three behaviours: a stationary state, undulatory swimming and steady spinning, where the swimming results from a self-oscillatory instability through a Hopf bifurcation. By conducting a linear stability analysis incorporating an elastohydrodynamic model, we theoretically predict the growth rates and critical conditions, which agree well with the numerical counterparts. We also propose a reduced model system consisting of a minimal elastic structure which reproduces the EEH instability. The elasto-viscous response of the composite structure is able to transform the pitchfork bifurcation into a Hopf bifurcation, leading to self-oscillation. Our results imply a new way of harnessing elastic media to engineer self-oscillations, and more generally, to manipulate and diversify the bifurcations and the corresponding instabilities. These ideas will be useful in designing soft, environmentally adaptive machines.

Tables

Table 1. Table 1: Parameters for the six cases chosen to validate numerical results against the theoretical force F z f → p subscript superscript 𝐹 → f p 𝑧 F^{\mathrm{f}\rightarrow\mathrm{p}}_{z} , torque Γ x f → p | V evaluated-at subscript superscript Γ → f p 𝑥 V \Gamma^{\mathrm{f}\rightarrow\mathrm{p}}_{x}|_{\mathrm{V}} and their leading-order counterparts. Bold entries indicate the difference with the reference, case 1 1 1 .

	$\tilde{\theta }$	$\delta$	${\sigma }_i$	$\overline{\mu }$
Case 1 (reference)	${10}^{-3}$	$0.3$	$0.2$	${10}^3$
Case 2	$0.1$	$0.3$	$0.2$	${10}^3$
Case 3	${10}^{-3}$	$0.8$	$0.2$	${10}^3$
Case 4	${10}^{-3}$	$0.3$	$\mathrm{𝟐}$	${10}^3$
Case 5	${10}^{-3}$	$0.3$	$0.2$	${\mathrm{𝟏𝟎}}^{\mathrm{𝟐}}$
Case 6	${10}^{-3}$	$0.3$	$0.2$	${\mathrm{𝟏𝟎}}^{\mathrm{𝟒}}$

Table 2. Table 2: Dimensionless critical electric field ℰ ¯ cri superscript ¯ ℰ cri \bar{\mathcal{E}}^{\mathrm{cri}} above which the composite object of a dielectric sphere of size ratio α 𝛼 \alpha and a rigid rod undergoes the electrohydrodynamic instability, the slenderness ϵ sl = 0.01 subscript italic-ϵ sl 0.01 \epsilon_{\mathrm{sl}}=0.01 .

$\alpha$	0.1	0.3	0.5	0.7	0.9
${\overline{ℰ}}^{\mathrm{cri}}$	5.278	1.803	1.348	1.202	1.136

Equations

$${\mathcal{P}}_i^{\infty }={\chi }^{\infty }{E}_i.$$

$$\frac{\partial {\mathcal{P}}_i}{\partial t}=-\frac{1}{{\tau }_{\mathrm{MW}}}\left[{\mathcal{P}}_i-\left({\chi }^0-{\chi }^{\infty }\right)E_i\right],$$

$$\kappa =\frac{R+2}{S+2}$$

$$E^{\mathrm{cri}}=\sqrt{\frac{2{\sigma }_{\mathrm{s}}\mu (R+2{)}^2}{3{ϵ}_{\mathrm{s}}^2(S-R)}}.$$

$${\mathbf{F}}^{\mathrm{f}\to \mathrm{p}}-{\beta }_{\mathrm{drag}}\mathbf{U}$$

$${\mathbf{\Gamma }}^{\mathrm{f}\to \mathrm{p}}+{\mathbf{\Gamma }}^{\mathrm{elec}}-{\gamma }_{\mathrm{drag}}\mathbf{\Omega }$$

$${\mathbf{\Gamma }}^{\mathrm{elec}}=\text{mathbfcal}{P}^{\mathrm{total}}\times \mathbf{E}=\text{mathbfcal}{P}^{\infty }\times \mathbf{E}+\text{mathbfcal}P\times \mathbf{E}=\text{mathbfcal}P\times \mathbf{E},$$

$$\frac{\partial \theta }{\partial t}$$

$$\frac{\partial \varphi }{\partial t}$$

$$\frac{\partial \psi }{\partial t}$$

$$\frac{\partial {\mathcal{P}}_N}{\partial t}+\frac{\partial \psi }{\partial t}{\mathcal{P}}_Q$$

$$\frac{\partial {\mathcal{P}}_Q}{\partial t}-\frac{\partial \psi }{\partial t}{\mathcal{P}}_N$$

$$\frac{\partial {\mathcal{P}}_3}{\partial t}$$

$$\frac{\partial \theta }{\partial \bar{t}}$$

$$\frac{\partial \varphi }{\partial \bar{t}}$$

$$\frac{\partial \psi }{\partial \bar{t}}$$

$$\bar{\eta }={\alpha }^3\bar{\mu },$$

$$\bar{\mu }=\frac{8\pi \mu L^4}{D{\tau }_{\mathrm{s}}}$$

$$\frac{\partial {\bar{\mathcal{P}}}_N}{\partial \bar{t}}$$

$$\frac{\partial {\bar{\mathcal{P}}}_Q}{\partial \bar{t}}$$

$$\frac{\partial {\bar{\mathcal{P}}}_3}{\partial \bar{t}}$$

$${\overline{\text{mathbfcal}P}}^{\infty }$$

$${\overline{\text{mathbfcal}P}}_{\mathrm{ini}}$$

$${\bar{\mathbf{F}}}^{\mathrm{f}\to \mathrm{p}}-3\alpha \bar{\mu }\bar{\mathbf{U}}/4$$

$${\bar{\mathbf{\Gamma }}}^{\mathrm{f}\to \mathrm{p}}+\overline{\text{mathbfcal}P}\times \bar{\mathbf{E}}-\bar{\eta }\bar{\mathbf{\Omega }}$$

$$8\pi \mu \left({\mathbf{r}}_t-{\mathbf{U}}^{\infty }\right)=c\left(\mathbf{I}+{\mathbf{r}}_s{\mathbf{r}}_s\right)\cdot \mathbf{f},$$

$$c=1+2\log {ϵ}_{\mathrm{sl}}<0.$$

$$\mathbf{f}(s)=-{\left(T(s){\mathbf{r}}_s\right)}_s+D{\mathbf{r}}_{ssss},$$

$$\bar{\mu }{\bar{\mathbf{r}}}_{\bar{t}}=-2c{\bar{T}}_{\bar{s}}{\bar{\mathbf{r}}}_{\bar{s}}-c\bar{T}{\bar{\mathbf{r}}}_{\bar{s}\bar{s}}+c{\bar{\mathbf{r}}}_{\bar{s}\bar{s}\bar{s}\bar{s}}+c\left({\bar{\mathbf{r}}}_{\bar{s}}\cdot {\bar{\mathbf{r}}}_{\bar{s}\bar{s}\bar{s}\bar{s}}\right){\bar{\mathbf{r}}}_{\bar{s}}.$$

$$2c{\bar{T}}_{\bar{s}\bar{s}}-c\bar{T}{\bar{\mathbf{r}}}_{\bar{s}\bar{s}}\cdot {\bar{\mathbf{r}}}_{\bar{s}\bar{s}}=-7c{\bar{\mathbf{r}}}_{\bar{s}\bar{s}}\cdot {\bar{\mathbf{r}}}_{\bar{s}\bar{s}\bar{s}\bar{s}}-6c{\bar{\mathbf{r}}}_{\bar{s}\bar{s}\bar{s}}\cdot {\bar{\mathbf{r}}}_{\bar{s}\bar{s}\bar{s}}-\bar{\mu }{\beta }_{\mathrm{p}}\left(1-{\bar{\mathbf{r}}}_{\bar{s}}\cdot {\bar{\mathbf{r}}}_{\bar{s}}\right),$$

$${\bar{\mathbf{r}}}_{\bar{s}\bar{s}}$$

$$\bar{T}$$

$$\bar{\mathbf{r}}$$

$${\bar{\mathbf{r}}}_{\bar{s}}$$

$$\frac{\mathrm{d}{\bar{\mathbf{x}}}_{\mathrm{p}}}{\mathrm{d}\bar{t}}$$

$$\frac{\mathrm{d}{\mathbf{e}}_{\mathrm{p}}}{\mathrm{d}\bar{t}}$$

$${\bar{\mathbf{F}}}^{\mathrm{f}\to \mathrm{p}}$$

Full text

March 17, 2024

1 Introduction

Active matter has been attracting much interest from a broad range of research communities (Ramaswamy, 2010; Cates & MacKintosh, 2011; Marchetti et al., 2013; Needleman & Dogic, 2017). At the micron scale, active matter consists of a large number of active agents that are able to convert energy to achieve directed or persistent motions, which include those of living microorganisms, synthetic micro-robots, biopolymers such as actin filaments, etc. The motions of these active agents are attributed to a wide range of mechanisms (Lauga & Powers, 2009; Marchetti et al., 2013; Alapan et al., 2019), e.g. one of the most common strategies adopted by natural and artificial micro-swimmers lies in the beating and wiggling of slender structures such as cilia and filaments, which are hair-like slender microscale structures that play an important role in various biological processes (Fawcett, 1961), such as swimming, pumping, mixing, cytoplasmic streaming, etc. The biological organelles deliver these functionalities by performing rhythmic, wave-like motions.

To achieve persistent motions, cyclic or oscillatory motions are needed, yet, the mechanism underlying the emergence of such oscillations remains unclear. Two major hypotheses, geometric feedback (Brokaw, 1971, 2009; Riedel-Kruse et al., 2007; Sartori et al., 2016; Hines & Blum, 1983; Hilfinger et al., 2009) and “flutter” or buckling instability (Bayly & Dutcher, 2016; De Canio et al., 2017; Ling et al., 2018; Hu & Bayly, 2018; Fatehiboroujeni et al., 2018), have been raised based on theory and/or simulations: the first hypothesis assumes that a time-dependent dynein activity (switching on/off or modulation) is necessary to trigger the oscillations; the second one suggests that a steady point force or force distributions acting along the axial direction of a flexible filament can trigger its oscillatory motion through a “flutter” or buckling instability. These forces are in fact called the “follower force” in the mechanics literature (Pflüger, 1950; Ziegler, 1952; Herrmann & Bungay, 1964). Because the follower force was initially invented theoretically and assumed to be always tangential to the slender structure regardless of its time-dependent deformation, it was demonstrated only mathematically and has been considered impractical (Koiter, 1996). However, it was recently realised experimentally on a metre-scale rod (Bigoni et al., 2018).

To drive the oscillations of artificial cilia and filaments of micron scale, different methods that exploit magnetic (Dreyfus et al., 2005; Singh et al., 2005; Evans et al., 2007; Livanovičs & Cēbers, 2012; Hanasoge et al., 2017; Huang et al., 2019) , electrostatic (den Toonder et al., 2008), piezoelectric (Kieseok et al., 2009), optical (van Oosten et al., 2009) and hydrogel-based actuations (Sidorenko et al., 2007; Masuda et al., 2013) have been developed. Nonetheless, these practises relied on a time-dependent power source, except for the self-oscillation of polymer brushes triggered by the Belousov-Zhabotinsky reaction (Masuda et al., 2013). This reaction-based beating shares with other biological processes, such as mammalian otoacoustic emissions (Gold, 1948; Kemp, 1979) and glycolysis (Sel’kov, 1968) the same feature: self-oscillation, that is generating and sustaining a periodic motion based on a power source without a corresponding periodicity (Jenkins, 2013).

In our recent work (Zhu & Stone, 2019), we proposed a chemical-reaction-free and follower-force-free strategy to engineer the self-oscillations of artificial structures by employing a time-independent, uniform electric field. We reported an elasto-electro-hydrodynamic (EEH) instability based on the Quincke rotation (QR) instability, and utilised it to drive various motions of an object composed of a dielectric spherical particle with an attached elastic filament. In this work, we will present in detail the setup and the mathematical description of the new EEH problem. First, we numerically solve the system coupling the electrohydrodynamics of the particle in a dielectric fluid and the elastohydrodynamics of the filament in a viscous fluid. We identify the emergence of the EEH instability that produces the self-oscillation of the composite object. The oscillations in turn cause the object to translate. Then, we perform a linear stability analysis (LSA) incorporating an elastohydrodynamic model to predict the onset of self-oscillatory instability. Finally, we propose a minimal model that can reproduce the similar EEH instability.

We describe the setup and governing equations of the EEH problem in § 2, and demonstrate the numerical results in § 3. The elastohydrodynamic model and LSA are shown in § 4, followed by § 5 illustrating the minimal model. Finally, we conclude our observations and provide some discussions in § 6.

2 Problem setup and mathematical formulations

We consider a weakly conducting dielectric spherical particle of radius $A$, which has attached an inextensible elastic filament of contour length $L$. The filament is cylindrical with a constant cross-section of radius $a$, and its slenderness is $\epsilon_{\mathrm{sl}}=a/L\ll 1$. We fix $\epsilon_{\mathrm{sl}}=0.01$ in this work. The composite object is subject to a time-independent uniform electric field $\mathbf{E}=E\mathbf{e}_{z}$ (see figure 1), where $\mathbf{e}_{z}$ is the $z$-direction basis vector of the laboratory coordinates system $\mathbf{e}_{xyz}$. The centreline of the filament is described by $\mathbf{r}\left(s,t\right)$, where $s$ indicates the arclength. The base J ($s=0$) of the filament is clamped at the particle surface, namely, the tangent vector $\partial\mathbf{r}/\partial s|_{s=0}=-\mathbf{e}_{\mathrm{p}}$ at the base always passes through the particle centre P, regardless of its deformation and the orientation vector $\mathbf{e}_{\mathrm{p}}$ of the particle. The size ratio between the particle and filament is $\alpha=A/L$. We consider only the bending deformation of the filament with a bending stiffness of $D=\pi a^{4}Y/4$, where $Y$ denotes Young’s modulus.

The composite object is immersed in a dielectric solvent fluid with dynamic viscosity $\mu$. The electrical conductivity and absolute permittivity of the solvent are $\sigma_{\mathrm{s}}$ and $\epsilon_{\mathrm{s}}$, respectively, and those of the particle are $\sigma_{\mathrm{p}}$ and $\epsilon_{\mathrm{p}}$; $R=\sigma_{\mathrm{p}}/\sigma_{\mathrm{s}}$ and $S=\epsilon_{\mathrm{p}}/\epsilon_{\mathrm{s}}$ indicate the ratios. The terms $\tau_{\mathrm{s}}=\epsilon_{\mathrm{s}}/\sigma_{\mathrm{s}}$ and $\tau_{\mathrm{p}}=\epsilon_{\mathrm{p}}/\sigma_{\mathrm{p}}$ denote the charge relaxation time of the solvent and particle, respectively. These electrical properties are important to the induced QR electrohydrodynamic instability that is critical to the dynamics in this paper. Their values are based on experiments (Brosseau et al., 2017), where $R=2.3\times 10^{-7}$ and $S=0.84$ are fixed in this work. Though the filament will also be polarised like the particle, the induced electric torque on the filament will be much weaker than that on the particle (see § 6 for a detailed discussion). We thus do not consider the electrohydrodynamics of the filament in this work.

2.1 Assumptions

The numerical simulations are carried out by invoking several assumptions. Motivated by biomimetic applications at the micron scale, we neglect the inertia of the fluid and particle. The fluid motion is therefore governed by the Stokes equations, and the particle satisfies instantaneous force- and torque-free conditions. The movement of the composite object is constrained to be planar, such that the particle centre P and filament position $\mathbf{r}(s,t)$ are in the $yz$-plane.

We adopt the local resistive-force theory (Batchelor, 1970) to calculate the hydrodynamic forces on the filament. We further ignore the hydrodynamic interactions between the particle and the filament. In the elastohydrodynamic model developed for LSA, we also assume that the filament undergoes weak deformation.

2.2 Electrohydrodynamics of the particle

When a dielectric particle in a dielectric solvent is exposed to an electric field, the interface of the particle will be polarised. The total induced dipole $\mathbfcal{P}^{\mathrm{total}}$ consists of an instantaneous part $\mathbfcal{P}^{\infty}$ and a retarding part $\mathbfcal{P}$, viz. $\mathbfcal{P}^{\mathrm{total}}=\mathbfcal{P}^{\infty}+\mathbfcal{P}$. Both vectors are defined by three components, $\mathcal{P}^{\infty}_{i}$ and $\mathcal{P}_{i}$ ($i=1...3$) in the reference frame $\mathbf{e}_{123}$ that rotates with the particle (see figure 2). For a homogeneous spherical particle, its Maxwell-Wagner polarisation time $\tau_{\mathrm{MW}}$, and low- and high-frequency susceptibilities, $\chi^{0}$ and $\chi^{\infty}$, respectively, are isotropic, hence the $i$-th component of the instantaneous dipole $\mathbfcal{P}^{\infty}$ is

[TABLE]

In the rotating reference frame of the particle, the retarding dipole $\mathbfcal{P}$ is governed by (Tsebers, 1980b; Cēbers et al., 2000)

[TABLE]

where

[TABLE]

and $\tau_{\mathrm{MW}}=\tau_{\mathrm{s}}/\kappa$. It is well known that when the charge relaxation time $\tau_{\mathrm{p}}$ of the particle is larger than that of the solvent $\tau_{\mathrm{s}}$, i.e., $R/S<1$, $\mathbfcal{P}$ is oriented opposite to the electric field. This directional misalignment is the necessary condition for the electro-rotation of the particle, the so-called Quincke rotation (Quincke, 1896), which occurs when, in addition, the strength $E$ of the electric field is above a critical value $E^{\mathrm{cri}}$ derived theoretically as (Jones, 1984; Brosseau et al., 2017) (see appendix B)

[TABLE]

We do not consider the hydrodynamic interactions between the spherical particle and filament, hence the dynamics of the particle can be obtained by using its translational and rotational mobility factors. Assuming that the particle rotates at angular velocity $\boldsymbol{\Omega}$ about its centre P, which translates at velocity $\mathbf{U}$, the force and torque balances on the particle give

[TABLE]

where $\mathbf{F}^{\mathrm{f}\rightarrow\mathrm{p}}$ denotes the force exerted by the filament on the particle, $\boldsymbol{\Gamma}^{\mathrm{f}\rightarrow\mathrm{p}}$ the torque with respect to the particle centre P, and $\beta_{\mathrm{drag}}=6\pi\mu A$ and $\gamma_{\mathrm{drag}}=8\pi\mu A^{3}$ are the translational and rotational drag coefficients of a sphere in the creeping flow, respectively. Also, $\boldsymbol{\Gamma}^{\mathrm{elec}}$ is the electric torque on the particle with respect to its centre P, that is

[TABLE]

where $\mathbfcal{P}^{\infty}\times\mathbf{E}\equiv\mathbf{0}$ for an isotropic sphere because $\mathcal{P}^{\infty}_{i}$ linearly scales with $E_{i}$ in each direction by the same factor $\chi^{\infty}$ (see equation (1)). It is worth noting that $\mathbfcal{P}^{\infty}\times\mathbf{E}\neq\mathbf{0}$ for ellipsoidal particles, where the factor $\chi^{\infty}$ is direction dependent (Cēbers et al., 2000; Brosseau et al., 2017). The translation of the particle is driven by the elastic force exerted by the filament, which is balanced by the viscous drag, while the rotational motion of the particle is determined by the balance between the elastic, electric and hydrodynamic torques.

The orientation of the particle $\mathbf{e}_{\mathrm{p}}$ is defined as the direction from the filament base J towards the particle centre P, where $\mathbf{e}_{3}$ of the particle-based reference frame coincides with $\mathbf{e}_{\mathrm{p}}$. We have found it convenient to use the proper Euler angles $[\theta,\phi,\psi]$, see figure 2. Here, $\mathbfcal{P}$ is decomposed into $\mathbfcal{P}=\mathcal{P}_{3}\mathbf{e}_{3}+\mathcal{P}_{N}\mathbf{N}+\mathcal{P}_{Q}\mathbf{Q}$, where $\mathbf{N}$ indicates the nodal line direction and $\mathbf{Q}=\mathbf{e}_{3}\times\mathbf{N}$. This decomposition applies to other vectorial variables such as $\mathbf{E}$. We constrain $\mathbfcal{P}$ onto the $yz$-plane, hence $\phi=\psi\equiv 0$ and $\theta$ is the only angle indicating the orientation $\mathbf{e}_{\mathrm{p}}$; additionally, $\mathcal{P}_{N}=0$ and $\mathbf{e}_{x}=\mathbf{N}$. For the sake of completeness, we first derive the governing equations for a general situation without these constraints.

Using the torque-free condition equation (5b), we obtain the governing equations for $[\theta,\phi,\psi]$,

[TABLE]

where $E_{3}=E\cos{\theta}$ and $E_{Q}=E\sin{\theta}$. The governing equations for $[\mathcal{P}_{N},\mathcal{P}_{Q},\mathcal{P}_{3}]$ are (Cēbers et al., 2000)

[TABLE]

We choose the charge relaxation time of the solvent $\tau_{\mathrm{s}}$ as the characteristic time, $L/\tau_{\mathrm{s}}$ the characteristic velocity, and $E^{\mathrm{cri}}$ and $D/(LE^{\mathrm{cri}})$ the characteristic strength of the electrical field and polarisation dipole, respectively. Using $\;\bar{}\;$ to indicate the dimensionless variables hereafter, the dimensionless equations for the Euler angles $[\theta,\phi,\psi]$ are

[TABLE]

as derived in Cēbers et al. (2000) in the absence of the elastic torque $\bar{\boldsymbol{\Gamma}}^{\mathrm{f}\rightarrow\mathrm{p}}$. Here,

[TABLE]

with

[TABLE]

defined as the elasto-electro-viscous (EEV) parameter. The dimensionless governing equations for $[\bar{\mathcal{P}}_{N},\bar{\mathcal{P}}_{Q},\bar{\mathcal{P}}_{3}]$ following from equation (8) are

[TABLE]

where $\kappa=(R+2)/(S+2)$ as defined in equation (3). We slightly perturb the instantaneous polarisation $\bar{\mathbfcal{P}}^{\infty}$ and use it as the initial value $\bar{\mathbfcal{P}}_{\mathrm{ini}}$ of $\bar{\mathbfcal{P}}$, where

[TABLE]

with $\epsilon_{\mathcal{P}}=\mathcal{O}(10^{-4})-\mathcal{O}(10^{-3})$.

The dimensionless force- and torque-free conditions are

[TABLE]

Since we constrain the motion of the composite object and the induced dipole $\bar{\mathbfcal{P}}$ to the $yz$-plane, we solve equations (9a), (12b) and (12c) for $\theta$, $\bar{\mathcal{P}}_{Q}$ and $\bar{\mathcal{P}}_{3}$, where the last term $\frac{\partial\psi}{\partial t}\bar{\mathcal{P}}_{N}$ in equation (12b) vanishes.

2.3 Elastohydrodynamics of the filament

We describe here the elastohydrodynamic equations for the filament. By employing the slender body theory (SBT) considering the leading-order local drag (Batchelor, 1970), the relation between the velocity $\mathbf{r}_{t}$ of the filament centreline and the force per unit length exerted by the fluid onto the filament $\mathbf{f}(s,t)$ is

[TABLE]

where $\mathbf{U}^{\infty}$ is the underlying flow velocity (background or imposed flow velocity) at $\mathbf{r}(s,t)$ and $\mathbf{U}^{\infty}=\mathbf{0}$ in this work; the subscripts $t$ and $s$ denote the partial derivatives with respect to $t$ and $s$, respectively and

[TABLE]

The filament is assumed to be described by the Euler–Bernoulli constitutive law, and because the elastic force balances the hydrodynamic force anywhere on the centreline, we obtain

[TABLE]

where $T(s,t)$ denotes the line tension, which acts as a Lagrangian multiplier to guarantee the inextensibility of the filament, i.e., $\mathbf{r}_{s}\cdot\mathbf{r}_{s}\equiv 1$.

By substituting equation (17) into equation (15), and choosing $L$ and $D/L^{2}$ as the characteristic length and force, respectively, we obtain the dimensionless equations for $\bar{\mathbf{r}}(\bar{s},\bar{t})$,

[TABLE]

The dimensionless equation for $\bar{T}\left(\bar{s}\right)$ reads,

[TABLE]

where the last term on the right-hand side $-\bar{\mu}\beta_{\mathrm{p}}\left(1-\bar{\mathbf{r}}_{\bar{s}}\cdot\bar{\mathbf{r}}_{\bar{s}}\right)$ is an extra (numerical) penalisation term introduced (Tornberg & Shelley, 2004; Li et al., 2013) to preserve the local inextensibility constraint $\bar{\mathbf{r}}_{\bar{s}}\cdot\bar{\mathbf{r}}_{\bar{s}}\equiv 1$; $\beta_{\mathrm{p}}=100$ is adopted in our simulations. The boundary conditions (BCs) for $\bar{\mathbf{r}}(\bar{s},\bar{t})$ and $\bar{T}(\bar{s},\bar{t})$ at the free end $\bar{s}=1$ are

[TABLE]

The BCs at the clamped end $\bar{s}=0$ couple the elastohydrodynamics and electrohydrodynamics, as will be described next.

2.4 Elasto-electro-hydrodynamic coupling

The electrohydrodynamics of the dielectric particle in a dielectric solvent and the elastohydrodynamics of the flexible filament in a viscous fluid are coupled via, first the BCs of $\bar{\mathbf{r}}(\bar{s},\bar{t})$ and $\bar{T}(\bar{s},\bar{t})$ at the filament base $\bar{s}=0$, and second the elastic force $\bar{\mathbf{F}}^{\mathrm{f}\rightarrow\mathrm{p}}(\bar{t})$ and torque $\bar{\boldsymbol{\Gamma}}^{\mathrm{f}\rightarrow\mathrm{p}}(\bar{t})$ exerted by the filament on the particle (equation (14).

The BCs at the filament base $\bar{s}=0$ are

[TABLE]

where $\bar{\mathbf{x}}_{\mathrm{p}}(\bar{t})$ denotes the dimensionless position of the particle centre P. Equations (21a) and  (21b) imply, respectively, that the filament base $\bar{s}=0$ is exactly on the particle surface, and the filament tangent vector at $\bar{s}=0$ always passes through the particle centre. Moreover, $\bar{\mathbf{x}}_{\mathrm{p}}(\bar{t})$ and $\mathbf{e}_{\mathrm{p}}(\bar{t})$ are connected to the particle kinematics through

[TABLE]

where $\bar{\mathbf{U}}(\bar{t})$ is linked to equation (14a) and $\bar{\boldsymbol{\Omega}}(\bar{t})$ to equation (9). The coupling is completed by the computation of $\bar{\mathbf{F}}^{\mathrm{f}\rightarrow\mathrm{p}}$ and $\bar{\boldsymbol{\Gamma}}^{\mathrm{f}\rightarrow\mathrm{p}}$,

[TABLE]

For completeness, we write the BC for the tension $\bar{T}$ at the filament base $\bar{s}=0$

[TABLE]

3 Numerical results

In the original QR phenomenon (without a filament), the particle rotates when the dimensionless electric field is above $1$, namely, $\bar{E}\geq 1$. We hereby investigate the influence of the bending stiffness of the filament by varying $\bar{\mu}$, where we fix the electric field $\bar{E}=1.5$ at which an individual particle undergoes steady QR. We fix the size ratio $\alpha=0.3$ in this section.

We observe that the composite object exhibits three $\bar{\mu}$-dependent scenarios, indicated by the time evolution of the rotational velocity $\bar{\Omega}$ shown in figure 3. When $\bar{\mu}=600$ (figure 3a), $\bar{\Omega}$ decays dramatically and eventually becomes zero, indicating that the object relaxes to a stationary state. Increasing $\bar{\mu}$ to $635$ (figure 3b), the time evolution of $\bar{\Omega}$ features two phases: in the initial phase (cyan domain), it grows rapidly due to self-oscillation; in the second phase (green domain), it reaches a time-periodic state with a constant amplitude of approximately $0.1$. The third type of response is illustrated by $\bar{\mu}=2000$, where $\bar{\Omega}$ eventually approaches a steady value around $-0.6$.

We further scrutinise the $\bar{\mu}=635$ case. The close-up views of the initially rapidly growing phase (cyan domain) and saturated time-periodic phase (green domain) are shown in figure 4a and b, respectively. The red curve connecting the local peaks $\bar{\Omega}^{\mathrm{lpk}}$ of $\bar{\Omega}$ implies an exponential growth of $\bar{\Omega}$ in time. This trend is confirmed by the linear relationship between $\log{\bar{\Omega}^{\mathrm{lpk}}}$ and $\bar{t}$ shown in the inset of figure 4a. The time-periodic phase enlarged in figure 4b reveals its sinusoidal-like variation characterised by fore-aft temporal symmetry. Six times within one period of this phase are marked, with their corresponding positions and orientations of the particle, and the profiles of the filament depicted in figure 4c. The oscillating particle drives the filament to wiggle, because the filament is clamped onto the particle. The wiggling filament provides thrust to the whole object, as a natural resemblance to a biological appendage. Consequently, the object achieves locomotion, following a wave-like trajectory (figure 4d). The wavy path is tightly packed near $\bar{t}=0$, implying the slow motion of the object undergoing small-amplitude oscillation in the initial phase.

We observe that when $\bar{\mu}$ lies in the self-oscillating regime, the time evolution of $\bar{\Omega}$ varies with $\bar{\mu}$. As shown in figure 5 for $\bar{\mu}=825$ and $1000$, for a larger $\bar{\mu}$ it takes fewer time periods for the perturbation to reach its time-periodic state. In addition, that state clearly breaks fore-aft symmetry with increasing $\bar{\mu}$.

We next investigate the critical $\bar{\mu}$ values that separate the three regimes corresponding to the stationary, undulating and steady spinning states. Figure 6 displays the rotational velocity magnitude $\bar{\Omega}^{\mathrm{mag}}$ versus $\bar{\mu}$ for $\bar{E}=1.2$ (a), $1.5$ (b) and $1.7$ (c). When $\bar{\mu}\leq\bar{\mu}_{1}^{\mathrm{cri}}$, $\bar{\Omega}^{\mathrm{mag}}=0$ represents the fixed-point solution; when $\bar{\mu}\geq\bar{\mu}_{1}^{\mathrm{cri}}$, the non-zero $\bar{\Omega}^{\mathrm{mag}}$ representing the constant spinning speed corresponds to the asymmetric fixed-point solution; when $\bar{\mu}\in\left(\bar{\mu}_{1}^{\mathrm{cri}},\bar{\mu}_{2}^{\mathrm{cri}}\right)$, $\bar{\Omega}^{\mathrm{mag}}$ indicates the magnitude of the oscillating $\bar{\Omega}$ when it reaches a time-periodic state. We plot $\left(\bar{\Omega}^{\mathrm{mag}}\right)^{2}$ as a function of $\bar{\mu}$ in close proximity to $\bar{\mu}_{1}^{\mathrm{cri}}$ in figure 6d-f. The linear dependence of $\left(\bar{\Omega}^{\mathrm{mag}}\right)^{2}$ on $\bar{\mu}$ implies that the instability occurs at $\bar{\mu}_{1}^{\mathrm{cri}}$ through a Hopf bifurcation from where a limit-cycle solution emerges. Moreover, the $\bar{\Omega}^{\mathrm{mag}}(\bar{\mu})$ profile also indicates the supercritical nature of the Hopf bifurcation. On the other hand, a sudden jump of $\bar{\Omega}^{\mathrm{mag}}$ at $\bar{\mu}_{2}^{\mathrm{cri}}$ signifies a secondary bifurcation where the limit cycle shrinks to a fixed point or vice versa.

Having demonstrated that the composite object is able to achieve propulsion by self-oscillatory undulation, we naturally examine its propulsive performance. Shown in figure 7a, when the undulating swimmer reaches its time-periodic state, its trajectory resembles a periodic wave propagating along a straight direction (dashed arrow). We thus define the effective translational velocity $\bar{\mathcal{U}}$ of the swimmer as the propagation speed of the wave, that is $\bar{\mathcal{U}}=\bar{\mathcal{D}}/\left(\bar{T}_{2}-\bar{T}_{1}\right)$. This effective velocity $\bar{\mathcal{U}}$ exhibits a clear non-monotonic variation in $\bar{\mu}$; it reaches its maximum value at an optimal EEV number $\bar{\mu}=\bar{\mu}^{\mathrm{opt}}\approx 825$ and becomes zero when $\bar{\mu}\rightarrow\bar{\mu}_{1}^{\mathrm{cri}}$ and $\bar{\mu}\rightarrow\bar{\mu}_{2}^{\mathrm{cri}}$. Such a non-monotonic trend is expected, since when $\bar{\mu}$ is outside the self-oscillating regime $[\bar{\mu}_{1}^{\mathrm{cri}},\bar{\mu}_{2}^{\mathrm{cri}}]$, the object either remains stationary or spins steadily, resulting in no net locomotion. It is also worth noting that $\bar{\mathcal{U}}$ exhibits wavy variation near $\bar{\mu}_{2}^{\mathrm{cri}}$. In this regime, the filament is so deflected and the hydrodynamic interactions between the particle and filament can be reasonably strong due to the decreasing distance between them. Since our simulations do not consider the hydrodynamic interactions, hence it is not self-consistent to interrogate the data in detail in this regime.

Finally, we show in figure 7c the dependence of the optimal swimming condition, $\bar{\mu}^{\mathrm{opt}}$ and $\bar{\mathcal{U}}^{\mathrm{opt}}$, on the electric field strength $\bar{E}$. The optimal EEV number $\bar{\mu}^{\mathrm{opt}}$ decreases with $\bar{E}$ monotonically; in contrast, the optimal velocity $\bar{\mathcal{U}}^{\mathrm{opt}}$ displays a non-monotonic variation in $\bar{E}$, reaching a maximum value of approximately $6\times 10^{-3}$ at $\bar{E}\approx 1.55-1.6$. This non-monotonic trend is not surprising. In fact, self-oscillation of the composite object only emerges when $1<\bar{E}<\bar{\mathcal{E}}^{\mathrm{cri}}$, where $\bar{\mathcal{E}}^{\mathrm{cri}}$ represents the critical electric field above which the particle jointed with a rigid rod ($\bar{\mu}\rightarrow 0$) of the same length and slenderness will undergo the QR instability. Hence, when $\bar{E}\geq\bar{\mathcal{E}}^{\mathrm{cri}}$, the composite object will spin steadily but not self-propel regardless of the filament rigidity. On the other hand, when $\bar{E}\leq 1$, the extra anchored filament will further stabilise the original QR particle, hence the composite object will be stationary. We further note that the optimal translational velocity $\approx 6\times 10^{-3}$ is in the range $\left(1,15\right)\times 10^{-3}$ of the dimensionless speed of a magnetically driven flexible artificial flagellum (Dreyfus et al., 2005).

By analogy to the results in figure 6a-c, we show in figure 8 $\bar{\Omega}^{\mathrm{mag}}$ versus $\bar{E}$ as the bifurcation parameter for three EEV numbers $\bar{\mu}=500$, $2000$ and $8000$. A similar bifurcation diagram is identified: increasing $\bar{E}$ from zero, the stationary fixed point solution transits to a limit-cycle solution through a supercritical Hopf bifurcation at $\bar{E}_{1}^{\mathrm{cri}}$ (green star); that solution then jumps to a second fixed point solution (steady spinning) via a secondary bifurcation at $\bar{E}_{2}^{\mathrm{cri}}$ (magenta star). The original QR instability emerges at $\bar{E}=1$ (hollow square) through a supercritical pitchfork bifurcation (Turcu, 1987; Peters et al., 2005; Das & Saintillan, 2013). The filament manages to transform that bifurcation for an individual particle into a corresponding Hopf bifurcation leading to self-oscillation. It is not surprising that by increasing $\bar{\mu}$, the variation of $\bar{\Omega}^{\mathrm{mag}}$ for the composite object tends to recover that of the original QR corresponding to $\bar{\mu}\rightarrow\infty$.

It is evident that the elastic torque $\bar{\boldsymbol{\Gamma}}^{\mathrm{f}\rightarrow\mathrm{p}}$ plays an important role in the torque balance. We examine the time evolution of the $x$-component of the torques, namely the elastic $\bar{\Gamma}_{x}^{\mathrm{f}\rightarrow\mathrm{p}}$, hydrodynamic $\bar{\Gamma}_{x}^{\mathrm{hydro}}$ and electric $\bar{\Gamma}_{x}^{\mathrm{hydro}}$ torques in figure 9 when $\bar{\mu}=635$ and $\bar{E}=1.5$. The sum of the torques $\bar{\Gamma}_{x}^{\mathrm{total}}=\bar{\Gamma}_{x}^{\mathrm{f}\rightarrow\mathrm{p}}+\bar{\Gamma}_{x}^{\mathrm{hydro}}+\bar{\Gamma}_{x}^{\mathrm{elec}}=0$ implies that the torque balance is well satisfied numerically. Similar to the evolution of the rotational velocity, the torques exhibit exponential growth in the initial phase before approaching a time-periodic state. The torque balance in this state is further scrutinised in figure 9b. Realising the negative relation between $\bar{\Gamma}^{\mathrm{hydro}}_{x}$ and $\bar{\Omega}$, we notice that $\bar{\Gamma}^{\mathrm{f}\rightarrow\mathrm{p}}_{x}$ and $\bar{\Omega}$ have the same sign in the two highlighted periods emphasising when the elastic $\bar{\Gamma}^{\mathrm{f}\rightarrow\mathrm{p}}_{x}$ and $\bar{\Gamma}^{\mathrm{hydro}}_{x}$ hydrodynamic torque contributions have opposite signs. The in-phase behaviour of $\bar{\Gamma}^{\mathrm{f}\rightarrow\mathrm{p}}_{x}$ and $\bar{\Omega}$ is a clear signature of negative damping, or positive feedback that triggers the linear instability of self-oscillation (Jenkins, 2013).

4 Linear stability analysis

4.1 Linearisation about the stationary equilibrium

state

We perform LSA about the stationary equilibrium state of the composite particle when the filament is undeformed. In this section, we drop the bars for all of the dimensionless unknown variables (those over dimensionless parameters remain), unless otherwise specified. We linearise the governing equations of the particle orientation $\theta$, and the dipole components $[\mathcal{P}_{Q},\mathcal{P}_{3}]$. By incorporating into the LSA a theoretical model of the elasto-viscous response of the filament, we do not linearise the equations for the filament position $\mathbf{r}(s)$ and tension $T(s)$ as conducted in Guglielmini et al. (2012).

The state variables $[\theta,\mathcal{P}_{Q},\mathcal{P}_{3}]$ are decomposed into a base (equilibrium) state $[\hat{\theta},\hat{\mathcal{P}}_{Q},\hat{\mathcal{P}}_{3}]$ and a perturbation state $[\theta^{\prime},\mathcal{P}^{\prime}_{Q},\mathcal{P}^{\prime}_{3}]$, which satisfy

[TABLE]

The perturbation-state variables $[\theta^{\prime},\mathcal{P}^{\prime}_{Q},\mathcal{P}^{\prime}_{3}]$ are assumed to be infinitesimal in LSA.

By substituting $\boldsymbol{\Gamma}^{\mathrm{f}\rightarrow\mathrm{p}}=\mathbf{0}$ and $\frac{\partial}{\partial t}=0$ into equations (9a), (12b) and 12c, we obtain the base-state dipoles

[TABLE]

By substituting equations (25) and (26) into equations (9a), (12b) and (12c), and assuming small $\theta^{\prime}$, we derive the governing equations for the perturbation-state variables $[\theta^{\prime},\mathcal{P}^{\prime}_{Q},\mathcal{P}^{\prime}_{3}]$,

[TABLE]

Adopting the normal-mode approach, we assume that the perturbations vary exponentially in time with a complex rate $\sigma=\sigma_{r}+\mathrm{i}\sigma_{i}$, so $\left[\mathcal{P}^{\prime}_{Q},\mathcal{P}^{\prime}_{3},\theta^{\prime}\right]=\left[\Phi,\Pi,\Theta\right]\exp{(\sigma t)}$. Consequently, equation (27) can be reformulated to

[TABLE]

We note that, for a vanishing elastic torque $\Gamma^{\mathrm{f}\rightarrow\mathrm{p}}_{N}=0$ (no attached filament), equation (28) is characterised by two roots $\sigma_{1}=0$ and $\sigma_{2}=\kappa\left(\bar{E}^{2}-1\right)$, which describe the original QR instability; the first root represents the stationary state and the second indicates that the dimensionless threshold electrical field (scaled by $E^{\mathrm{cri}}$) required to trigger instability is $\bar{E}=1$. Note that $E^{\mathrm{cri}}$ in equation (4) is originally derived by balancing the electric and hydrodynamic torque (Jones, 1984) instead of conducting LSA (see appendix B for details). The two predictions exactly agree with each other.

4.2 Elastohydrodynamic model

Since the elastohydrodynamic equations are not linearised, we thus derive a theoretical expression for $\Gamma^{\mathrm{f}\rightarrow\mathrm{p}}_{N}(t)$ for the dispersion relation, equation (28).

We find $\Gamma^{\mathrm{f}\rightarrow\mathrm{p}}_{N}(t)$ by solving a separate elastohydrodynamic problem of the composite object undergoing a prescribed rotational oscillation characterised by $\theta^{\prime}=\tilde{\theta}(t)\exp{\left(\mathrm{i}\sigma_{i}t\right)}$, where $\tilde{\theta}(t)=\Theta\exp{\left(\sigma_{r}t\right)}$ indicates the angular oscillation amplitude. We do not consider the object’s translation near the onset of instability since any translation is negligible due to the small-amplitude oscillation. To simplify the algebra in the next steps, we set $\hat{\theta}=\pi/2$ without loss of generality as shown in figure 10, where the rest configuration (dashed curves) corresponds to when the particle centre P coincides with the origin O and the undeformed filament is aligned in the $\mathbf{e}_{y}$ direction. The rotational oscillation is executed about a pivot V that lies away from the origin by a dimensional distance $b$ on the $y$-axis, where $\beta=b/L$; the dimensional distance between V and J is $d=A-b$, so similarly

[TABLE]

The particle centre P (resp. filament base J) follows a trajectory of a circular arc that is centred at V and of radius $b$ (resp. $d$); both trajectories are symmetric about the $y$-axis. Note that $\beta$ is an unknown that is to be determined.

Near the onset of instability, the amplitude $\tilde{\theta}(t)$ varies much more slowly than the oscillation of $\theta^{\prime}$, viz. $\sigma_{r}\ll\sigma_{i}$. This allows us to assume that the amplitude $\tilde{\theta}=\Theta\exp{(\sigma_{r}t)}$ is quasi-steady, namely, $\theta^{\prime}$ at a particular time $t_{0}$ can be approximated by

[TABLE]

as an instantaneous configuration of a periodic signal with a prescribed amplitude $\Theta\exp{\left(\sigma_{r}t_{0}\right)}$ and frequency $\sigma_{i}$. This setup resembles the theoretical framework developed to address the so-called elastohydrodynamic problem II (Wiggins & Goldstein, 1998; Wiggins et al., 1998) of a filament with one of its ends undergoing straight, oscillatory translation. We adapt that framework for our configuration, whereas the filament end oscillates on a circular arc instead of on a straight path, as shown in figure 10b. Because the filament undergoes small-amplitude deformation, $|z_{y}|\ll 1$ and its tangent vector $\mathbf{r}_{s}\approxeq\mathbf{e}_{y}$. We also assume $T(s)\equiv 0$. The position $\mathbf{r}(t,s)$ of the filament centreline is $\mathbf{r}(s)=(\alpha+s)\mathbf{e}_{y}+z(t,s)\mathbf{e}_{z}$. The horizontal displacement of the filament base is of order $\mathcal{O}(\tilde{\theta}^{2})$ and can be neglected because $|\theta^{\prime}|\leq|\tilde{\theta}|\ll 1$. The base’s vertical oscillation is prescribed as

[TABLE]

where $\delta\tilde{\theta}$ represents the oscillation amplitude. Following Wiggins & Goldstein (1998) and Wiggins et al. (1998), the deflection of the filament is expressed by

[TABLE]

where

[TABLE]

and $h$ is a sum of four solutions

[TABLE]

with

[TABLE]

The four coefficients $c_{i}$ need to be determined by the BCs at the filament ends. In contrast to Wiggins & Goldstein (1998) and Wiggins et al. (1998) treating $z(s)$ as a real variable, we consider a complex $z(s)$. This allows us to obtain the complex torque consistent with the complex nature of the torque balance, equation (28).

The BCs for $h(s)$ at the free end $s=1$ are $h_{ss}=h_{sss}=0$. At the clamped end $s=0$, $h=1$ as a Dirichlet BC corresponding to the prescribed displacement; the other BC is more subtle. Because the filament orientation is orthogonal to the circular arc (see figure 10b), we have

[TABLE]

By substituting equation (32) into equation (36), we obtain the BC

[TABLE]

where $\delta$ is defined in equation (29). Knowing all the BCs of $h(s)$, we compute the four coefficients

[TABLE]

where $\Lambda=\xi^{2\mathrm{i}}+4\xi^{1+\mathrm{i}}+\xi^{2+2\mathrm{i}}+\xi^{2}+1$. Considering the small-amplitude deformation, the total force $\mathbf{F}$ exerted by the filament on the clamped end is along the vertical $\mathbf{e}_{z}$ direction. The torque $\boldsymbol{\Gamma}^{\mathrm{f}\rightarrow\mathrm{p}}|_{V}$ with respect to the pivot V and $\boldsymbol{\Gamma}^{\mathrm{f}\rightarrow\mathrm{p}}$ with respect to the particle centre P are along the $\mathbf{e}_{x}$ direction, so that the corresponding components of the force and torques are

[TABLE]

where

[TABLE]

Now, let us examine the denominator, $\Lambda$, of equation (39) whose five terms are in the form of $\xi^{q_{k}}$ ($k=1...5$), where $[q_{1},q_{2},q_{3},q_{4},q_{5}]=[2\mathrm{i},1+\mathrm{i},2+2\mathrm{i},2,0]$. Using equation (35), we express $\xi^{q_{k}}$ as

[TABLE]

where $\zeta_{k}=\exp{(z_{0}q_{k})}$ are

[TABLE]

We observe that the third term $\zeta_{3}^{\mathcal{L}}$ is larger than the rest in magnitude when $\mathcal{L}\geq 1$, dominating the second largest term by one order when $\mathcal{L}\geq 3$. Let us assume $\mathcal{L}\geq 3$ a priori, so that we can then approximate $\Lambda$ by $\zeta_{3}^{\mathcal{L}}$ in equation (39). By further extracting the leading-order terms of $\Lambda_{1}/\Lambda$, $\Lambda_{2}/\Lambda$ and $\Lambda_{3}/\Lambda$, we attain a simplified, leading-order expression for the force and torque (denoted by $\;\tilde{\;}\;$)

[TABLE]

The theoretical force $F^{\mathrm{f}\rightarrow\mathrm{p}}_{z}$, torque $\Gamma^{\mathrm{f}\rightarrow\mathrm{p}}_{x}|_{\mathrm{V}}$ and their leading-order counterparts $\tilde{F}^{\mathrm{f}\rightarrow\mathrm{p}}_{z}$ and $\tilde{\Gamma}^{\mathrm{f}\rightarrow\mathrm{p}}_{x}|_{\mathrm{V}}$ are validated against the numerical results for six cases spanning a wide range of parameters relevant to our study (see table. 1), where case 1 is the reference case and the other five vary a single parameter compared to case 1. Because the numerical force and torque are real quantities, the real parts of $F^{\mathrm{f}\rightarrow\mathrm{p}}_{z}$ (dashed curve) given by equation (39a), and its leading-order approximation $\tilde{F}^{\mathrm{f}\rightarrow\mathrm{p}}_{z}$ (dot-dashed curve) by equation (43a), are compared with the numerical data (solid curve) in figure 11. A similar comparison between the torques $\Gamma^{\mathrm{f}\rightarrow\mathrm{p}}_{x}|_{\mathrm{V}}$ and $\tilde{\Gamma}^{\mathrm{f}\rightarrow\mathrm{p}}_{x}|_{\mathrm{V}}$ is shown in figure 12.

We observe that the force $F^{\mathrm{f}\rightarrow\mathrm{p}}_{z}$ and torque $\Gamma^{\mathrm{f}\rightarrow\mathrm{p}}_{x}|_{\mathrm{V}}$ and their leading-order values agree with the numerical results quantitatively in all the cases except for case $5$, where the leading-order results deviate a little from the full expression and numerical results. This disagreement results from the violation of the assumption $\mathcal{L}\geq 3$ used to derive the leading-order expression, where $\mathcal{L}\approx 1.25$ for case $5$. This also implies that the leading-order predictions become less accurate at small $\bar{\mu}$ values.

For the validation purpose, $\delta=\alpha-\beta$ can be prescribed. However, for the model, $\delta$ needs to be determined using the force-free condition on the particle. The particle follows a circular arc on the other side of the pivot V, the $z$-component of the hydrodynamic force on the particle approximated by the Stokes’s law is

[TABLE]

Substituting equation (43a) and (44) into $\tilde{F}^{\mathrm{f}\rightarrow\mathrm{p}}_{z}+F_{z}^{\mathrm{f}\rightarrow\mathrm{p}}=0$, we obtain

[TABLE]

Using the leading-order torque $\tilde{\Gamma}^{\mathrm{f}\rightarrow\mathrm{p}}_{x}$ equation (43c), as the left-hand side torque of equation (28) (note that the nodal line direction $\mathbf{N}=\mathbf{e}_{x}$ when the orientation $\mathbf{e}_{\mathrm{p}}$ is restricted to the $yz$-plane), we obtain the governing equation for the transformed growth rate $\hat{\sigma}=\bar{\mu}\sigma$,

[TABLE]

where

[TABLE]

where $\log{\xi}$ can be written as

[TABLE]

4.3 Complex growth rates and onset of

instability

We solve equation (46) to obtain the transformed growth rate $\hat{\sigma}$, which facilitates a theoretical prediction of the onset of self-oscillatory instability. The growth rate $\hat{\sigma}=\hat{\sigma}_{r}+\mathrm{i}\hat{\sigma}_{i}$ depends on $\alpha$, $\epsilon_{\mathrm{sl}}$, $\kappa$, $\bar{\mu}$ and $E$, where we have fixed $\epsilon_{\mathrm{sl}}$ and $\kappa$. By writing $\hat{\sigma}_{i}=W^{4}$ and substituting it into equation (33), we obtain $\mathcal{L}=W/(-1-2\log{\epsilon_{\mathrm{sl}}})^{1/4}$. Here, $\mathcal{L}$ is a positive real number, so is $W$. By substituting equations (47) and (48) into equation (46), we derive a system of two-dimensional, nonlinear polynomial equations for $\hat{\sigma}_{r}$ and $W$ (see appendix C) and obtain its roots by employing the python driver phcpy (Verschelde, 2013; Otto et al., 2019) of a general-purpose solver PHCpack (Verschelde, 1997) for polynomial systems. Because $\hat{\sigma}=\bar{\mu}\sigma$, we obtain the real part $\sigma_{r}=\hat{\sigma}_{r}/\bar{\mu}$ and imaginary part $\sigma_{i}=W^{4}/\bar{\mu}$ of the complex growth rate $\sigma$.

We show $\sigma_{r}$ and $\sigma_{i}$ as a function of $\bar{\mu}$ in figure 13 for two electric fields $\bar{E}=1.15$ (a) and $\bar{E}=1.5$ (b), where $\alpha=0.3$. In both cases, the imaginary part $\sigma_{i}\left(\bar{\mu}\right)>0$ implying that the perturbation always decays/grows in an oscillatory manner. In contrast, the real part $\sigma_{r}$ increases with $\bar{\mu}$ monotonically from negative to positive values, indicating the critical condition $\sigma_{r}\left(\bar{\mu}_{1}^{\mathrm{cri}}\right)=0$ of the self-oscillatory instability. When $\bar{\mu}$ is smaller/larger than $\bar{\mu}_{1}^{\mathrm{cri}}$, the perturbation exhibits oscillatory decaying and growth. The LSA prediction of $\left(\sigma_{r},\sigma_{i}\right)$ agrees quantitatively with the numerical counterpart for the $\bar{E}=1.15$ case, and qualitatively for the $\bar{E}=1.5$ case.

We adopt a bi-section method to determine $\bar{\mu}_{1}^{\mathrm{cri}}$ as a function of $\left(\bar{E},\alpha\right)$, as shown in figure 14. For all $\alpha$ values, $\bar{\mu}_{1}^{\mathrm{cri}}$ decreases monotonically with $\bar{E}$. The theoretical and numerical predictions agree well with each other, especially in the high $\bar{\mu}$ regime. The agreement degenerates with decreasing $\bar{\mu}$. We infer that $\mathcal{L}$ becomes smaller when $\bar{\mu}$ decreases, hence this disagreement is mostly attributed to violating the $\mathcal{L}\geq 3$ assumption of the leading-order force/torque model for the LSA.

5 A minimal model to reproduce the EEH instability and

self-oscillation

To better unravel the physics underlying the EEH instability, we seek a minimal model reproducing this instability and the corresponding self-oscillation. By analogy to the multi-linker models (De Canio et al., 2017; Ling et al., 2018), we replace the elastic filament by two rigid cylindrical rods numbered $\#1$ and $\#2$ of equal length $\ell=L/2$ and equal radius $a$ of their cross sections, which are linked at $\mathrm{J}_{1}$ by a torsional spring with an elastic module of $K$ (see figure. 15). Rod $\#1$ is clamped at the sphere surface J, namely it always passes through the particle centre P, hence the displacement vectors $\overrightarrow{\mathrm{PJ}}$ and $\overrightarrow{\mathrm{P}\mathrm{J}_{1}}$ are opposite to the particle orientation $\mathbf{e}_{\mathrm{p}}$. Rod $\#2$ is oriented with respect to $\overrightarrow{\mathrm{P}\mathrm{J}_{1}}$ by an angle $\theta_{1}$, which is zero when the composite system is at rest.

Similar to the original setup, we assume that the motion of particle and the rods are restricted to the $yz$-plane. Further, no hydrodynamic interactions between the particle and rods, or between the rods are considered. The system consists of six unknowns: the translational velocity components $U_{y}(t)$ and $U_{z}(t)$ of the particle, the rotational velocity component $\frac{\mathrm{d}\theta(t)}{\mathrm{d}t}=\Omega(t)$ of the particle and $\frac{\mathrm{d}\theta_{1}(t)}{\mathrm{d}t}=\Omega_{1}(t)$ of rod $\#2$ with respect to rod $\#1$, and the polarisation vector components $\mathcal{P}_{Q}(t)$ and $\mathcal{P}_{3}(t)$. It is worth noting that compared to the classical QR particle, this minimal configuration only incorporates one extra degree of freedom, $\theta_{1}$, which indicates the deformation magnitude of the torsional spring.

We first derive the hydrodynamic force exerted on rod $\#1$ as

[TABLE]

and the torque about the particle centre P

[TABLE]

Likewise, the hydrodynamic force exerted on rod $\#2$ is

[TABLE]

and the hydrodynamic torque on rod $\#2$ about $\mathrm{J}_{1}$ is

[TABLE]

The torque-free condition on rod $\#2$ reads

[TABLE]

where $\mathbf{M}_{2}=-K\theta_{1}\mathbf{e}_{x}$ is the elastic moment exerted on rod $\#2$ by the torsional spring. The torque balance on the whole composite system about the particle centre P is

[TABLE]

We also need to impose the force-free condition on the whole composite object

[TABLE]

To close the system, we solve the governing equations (8b) and (8c) for $\mathcal{P}_{Q}$ and $\mathcal{P}_{3}$, where the second term $-\frac{\partial\psi}{\partial t}\mathcal{P}_{N}$ in equation (8b) disappears. We note that equations (53) and (54) indeed reflect the subtle interplay between the elastic, electric and hydrodynamic torques, which lead to the EEH instability-induced self-oscillation.

5.1 Nondimensionalization of the minimal model

We use the same characteristic scales as the original particle-filament configuration (see § 2) to nondimensionalise equations (53), (54) and (55), except that we substitute $D$ by $KL$, resulting in a slightly modified EEV parameter

[TABLE]

to be distinguished from $\bar{\mu}$ defined by equation (11) for the original setup. The dimensionless governing equations for $\bar{U}_{y}(\bar{t})$, $\bar{U}_{z}(\bar{t})$, $\theta(\bar{t})$ and $\theta_{1}(\bar{t})$ are

[TABLE]

where $c=1+2\log\epsilon_{\mathrm{sl}}$ and $\bar{\eta}=\alpha^{3}\breve{\mu}$ as given by equations (16) and (10), respectively.

The dimensionless equations for $\bar{\mathcal{P}}_{Q}$ and $\bar{\mathcal{P}}_{3}$ are

[TABLE]

with their initial values at $\bar{t}=0$

[TABLE]

5.2 Numerical and theoretical (LSA) results of the minimal

model

We solve equations (57) and (58) numerically using the MATLAB solver ‘ode15s’ for ordinary differential equations. Fixing the electric field $\bar{E}=1.5$ and size ratio $\alpha=0.3$, we show in figure 16a the $\breve{\mu}$-dependent magnitudes $\bar{\Omega}^{\mathrm{mag}}$ and $\bar{\Omega}^{\mathrm{mag}}_{1}$ of the rotational velocities of the particle and rod $\#2$, respectively, when the minimal composite object reaches its equilibrium configuration. This simple model reproduces the three characteristic behaviours of the original particle-filament system: stationary ($\breve{\mu}<\breve{\mu}_{1}^{\mathrm{cri}}\approx 2513$), wiggling ($\breve{\mu}_{1}^{\mathrm{cri}}<\breve{\mu}<\breve{\mu}_{2}^{\mathrm{cri}}\approx 4300$) and steady spinning ($\breve{\mu}>\breve{\mu}_{2}^{\mathrm{cri}}$). In the spinning state, $\bar{\Omega}^{\mathrm{mag}}_{1}=|\mathrm{d}\theta_{1}/\mathrm{d}\bar{t}|=0$ reflects a time-independent angle $\theta_{1}$ between the two rods, which adopt a steady “deformed” configuration representing a minimal model of the deformed filament.

Conducting an LSA for this minimal model, we find the closed-form expression of the complex growth rate $\sigma=\sigma_{r}+\mathrm{i}\sigma_{i}$ (see appendix D for details). The theoretical values of $\sigma_{r,i}$ versus $\breve{\mu}$ in case of $\left(\bar{E},\alpha\right)=\left(1.5,0.3\right)$ are depicted in figure 16b, as well as their numerical counterparts in the near $\breve{\mu}_{1}^{\mathrm{cri}}$ regime. The theoretical and numerical values of both $\sigma_{r}$ and $\sigma_{i}$ almost lie on top of each other, consequently, their predictions of $\breve{\mu}_{1}^{\mathrm{cri}}$ (when $\sigma_{r}=0$) agree. This superior agreement to the particle-filament system (figure 13) is expected, because the minimal model does not require an approximate model (see §4.2) for the elastic torque as the original case.

The LSA also indicates the emergence of another subtle critical EEV number $\breve{\mu}_{0}^{\mathrm{cri}}\approx 345$ (black cross in figure 16b): when $\breve{\mu}<\breve{\mu}_{0}^{\mathrm{cri}}$, the real part $\sigma_{r}$ of the growth rate is negative, accompanying a zero imaginary part, thus the perturbations diminish to zero monotonically; when $\breve{\mu}_{0}^{\mathrm{cri}}<\breve{\mu}<\breve{\mu}_{1}^{\mathrm{cri}}$, $\sigma_{r}<0$ but $\sigma_{i}>0$, the perturbations also die out but in a an oscillatory fashion. The former case corresponds to the non-negative quantity $\Sigma$ inside the square-root operator in equation (85) that naturally yields real solutions for $\sigma$ only. A similar structure of the solutions of $\sigma$ was reported in De Canio et al. (2017). Since the current work mainly addresses the EEH instability-induced self-oscillation, we do not pursue a detailed investigation in this stable, stationary regime.

6 Conclusions and discussions

Standard biomimetic practises commonly rely on an oscillating magnetic or electric field to produce the oscillatory motion of slender artificial structures. In contrast, we propose a strategy to achieve self-oscillation of artificial structures based on a time-independent, uniform electric field. By formulating and numerically solving an elasto-electro-hydrodynamic problem, this concept is illustrated by oscillating a composite object consisting of a weakly conducting dielectric spherical particle and an elastic filament immersed in a dielectric solvent.

Our strategy is grounded in the QR electrohydrodynamic instability phenomenon indicating that a weakly conducting dielectric particle suspended in a dielectric liquid of higher conductivity can undergo spontaneous rotation under a sufficiently strong DC electric field. For an individual spherical particle, this instability emerges through a supercritical pitchfork bifurcation resulting in steady rotation (Jones, 1984). By incorporating an elastic filament, we transform the pitchfork bifurcation into a Hopf bifurcation through which a self-oscillatory instability occurs (Zhu & Stone, 2019). This transformation is attributed to the elasto-viscous response of the filament providing an elastic torque to balance the electric and hydrodynamic torques. The elastic torque is in phase with the rotational velocity of the particle at certain time periods (see figure 9b). This in-phase behaviour results in negative damping (or positive feedback), hence leading to the onset of linear instability (Jenkins, 2013). We comment that such a transition from pitchfork to Hopf bifurcation was also identified by Tsebers (1980a) who observed oscillatory QR of ellipsoidal particles attributed to their anisotropic electric properties. It is also worth mentioning that the QR instability was utilised to study suspensions of artificial swimmers made of QR particles that achieved locomotion by rolling near a rigid solid boundary (Bricard et al., 2013). In addition, the recent work of Das & Lauga (2019) shows theoretically and numerically that a dielectric particle with particular geometrical asymmetry (e.g. a helix) under a DC electric field is able to convert QR into spontaneous translation in an unbounded domain.

We next recall the original experiments conducted by Quincke (1896), where the particle was hung by a silk thread and hence the particle rotated in the direction along the orientation of the thread. Quincke also noted an oscillatory behaviour as translated by Jones (1984)

“Quincke, with his spheres tethered to silk threads, had been forced to contend with periodic rotation, first in one direction and then in the other as the silk thread wound and unwound”.

We think that the “wound and unwound” motion manifested the self-oscillatory phenomenon, which is attributed to the torsional deformation of the silk thread. We speculate that Quincke probably regarded this observation as an experimental nuisance, thus did not pay attention to it nor did other researchers, except for one little-known preprint (Zaks & Shliomis, 2014) that recognised and modelled this torsional oscillation by considering a QR particle hung by a thread with torsional elasticity.

In this paper, we consider only the bending stiffness of the grafted filament and the whole composite object is freely suspended in the solvent. By applying an electric field stronger than the critical value corresponding to the onset of original QR instability, the composite object exhibits three distinct behaviours depending on the EEV number $\bar{\mu}$ (inversely proportional to the bending stiffness). When $\bar{\mu}\leq\bar{\mu}_{1}^{\mathrm{cri}}$, the object remains stationary, corresponding to a fixed-point solution; when $\bar{\mu}\geq\bar{\mu}_{2}^{\mathrm{cri}}$, the particle spins steadily towing a deformed filament, corresponding to an asymmetric fixed-point solution; when $\bar{\mu}\in\left(\bar{\mu}_{1}^{\mathrm{cri}},\bar{\mu}_{2}^{\mathrm{cri}}\right)$, the particle oscillates and the filament wiggles, leading the object to an undulatory locomotion. More specifically, instability occurs at $\bar{\mu}_{1}^{\mathrm{cri}}$ through a supercritical Hopf bifurcation, where the self-oscillatory motion represents a limit-cycle solution; at $\bar{\mu}_{2}^{\mathrm{cri}}$, a secondary bifurcation appears, and the oscillatory, limit-cycle solution jumps to the steadily spinning, fixed-point solution. By fixing the EEV number $\bar{\mu}$, bifurcation diagrams considering the electric field strength $\bar{E}$ as the control parameter revealed the same three scenarios (see figure 8).

We have also examined the propulsive performance of the micro object in the self-oscillating regime $\bar{\mu}\in(\bar{\mu}_{1}^{\mathrm{cri}},\bar{\mu}_{2}^{\mathrm{cri}})$. The trajectory of the object resembles a wave propagating along a straight path. The translational velocity of the object along this path varies in $\bar{\mu}$ non-monotonically (see figure 7c).

Motivated by the exponential temporal growth of the rotational velocity, we performed a LSA to predict theoretically the onset of the self-oscillatory instability. We have developed an elastohydrodynamic model to account for the elastic force and torque exerted by the filament on the particle, which closely matched the numerical counterparts. Incorporating this model into a standard LSA for the original QR particle, we derived the dispersion relationship of the new EEH problem. We thus calculated the complex growth rate $\sigma=\sigma_{r}+\mathrm{i}\sigma_{i}$ and identified the critical EEV number $\bar{\mu}_{1}^{\mathrm{cri}}$. Theoretical predictions of $\sigma$ (figure 13) and $\bar{\mu}_{1}^{\mathrm{cri}}$ (figure 14) agree well the numerical results, especially in the large $\bar{\mu}$ regime. However, the agreement becomes less satisfactory when $\bar{\mu}$ decreases because of the violation of an assumption used in the elastohydrodynamic model.

To unravel the EEH instability mechanism, we studied a minimal model system characterised by two rigid rods linked by a torsional spring to mimic the original filament. This substitution reduces the elastic element’s number of degrees of freedom to one. Numerical and LSA results demonstrated that the minimal model could exhibit the three elasticity-dependent behaviours: stationary, wiggling and steady spinning.

Following the comments of an anonymous referee, we hereby emphasise the difference between our work and other seemingly similar studies (Manghi et al., 2006; Qian et al., 2008; Coq et al., 2008), where a flexible slender structure (filament or rod) rotated in a viscous fluid and produced thrust because one of its ends was clamped to a constantly rotating base or actuated by a constant torque. This rotation results from forced oscillation characterised by a close correlation between the frequency of the power source and that of the resulting periodic motion. This forced-oscillatory periodic motion distinguishes itself from the self-oscillatory motion we observe, where the time-independent electric field as the power source lacks a frequency corresponding to that of the periodic motion.

The current work constrained the kinematics and electric polarisation vector of the particle to a plane in order to show a clean physical picture of the new EEH instability we identified. By removing these constraints, we anticipate the appearance of more complex and diverse three-dimensional behaviours featured by bi/multi-stability, hysteresis and even chaos (ellipsoidal particles were observed to exhibit chaotic QR (Tsebers, 1991)). We will report the results of the ongoing work in a future paper.

It is also worth mentioning the assumption of neglecting electrohydrodynamic effect of the filament. The electric torque exerted on a slender QR structure scales with $a^{2}L$ (Das & Lauga, 2019), and that on a sphere scales with $A^{3}$ (see equation (69)). By assuming that the filament and particle have similar dielectric properties and realising $\alpha=A/L=O(1)$, the ratio of the former to the latter torque is of the order of $\epsilon_{\mathrm{sl}}^{2}$. This comparison thus justifies the assumption, which also implies that no special attention needs to be paid in this context for the experimental realisation.

In conclusion, incorporating an elastic element to manipulate the electrohydrodynamic instability, we report an elasto-electro-hydrodynamic instability and use it for engineering self-oscillation of artificial structures. We anticipate that this idea of harnessing elastic media to control and diversify the bifurcation and the corresponding instability behaviour can be generalised to other stability phenomena and systems. As a result, different emerging instability behaviours can be utilised for diverse functionalities. This concept might inspire new approaches to design soft, reconfigurable machines that can morph and adapt to the environment.

Declaration of Interests. The authors report no conflict of interest.

Appendix A Numerical methods for the EEH problem

In this section, we describe the numerical methods to solve the EEH problem. For notation brevity, we remove all the bars over the unknown dimensionless variables henceforth. Following Tornberg & Shelley (2004), we use a finite difference scheme to discretise the filament centreline by a uniform grid of $N$ points. A typical value of $N=201$ is used in the simulations. Therefore, we have $N$ unknowns for the tension $T(s)$ and $3N$ for the coordinates $\mathbf{r}(s)=x\mathbf{e}_{x}+y\mathbf{e}_{y}+z\mathbf{e}_{z}$. The implementation considers a general three-dimensional motion of the filament, hence we also solve for $x(s)$ even though the motion is restricted to the $yz$-plane.

In contrast, for the restricted motion of the particle, we express its translational velocity $\mathbf{U}$ and rotational velocity $\boldsymbol{\Omega}$ as

[TABLE]

which yields three unknowns $U_{y}$, $U_{z}$ and $\Omega=\theta_{t}$ for the particle. Another two unknowns are $\mathcal{P}_{Q}$ and $\mathcal{P}_{3}$ governed by equations (12b) and (12c), respectively.

At the $(k+1)$-th time step, $T^{k+1}(s)$ is solved based on $\mathbf{r}^{k}(s)$ at the $k$-th time step. We then solve $[\mathbf{r}^{k+1}(s),U^{k+1}_{y},U^{k+1}_{z},\Omega^{k+1},\mathcal{P}^{k+1}_{Q},\mathcal{P}^{k+1}_{3}]$ in a coupled way, which consists of $3N+5$ unknowns. We adopt this coupled strategy to accurately preserve the clamped BC of the filament base J ($s=0$): first, the filament base is on the particle surface; second, the tangent vector $\mathbf{r}_{s}|_{s=0}$ at the base always passes through the particle centre P. This clamped BC of the filament is different from other configurations (Guglielmini et al., 2012; De Canio et al., 2017) where the base is stationary. In our case, the filament base is right on the particle surface, and is able to translate with and rotate about the particle centre.

We use the backward Euler formulation to approximate the particle centre $\mathbf{x}_{\mathrm{p}}^{k+1}$,

[TABLE]

where $\Delta t$ is the time step. The Dirichlet BC for $\mathbf{r}|_{s=0}$, namely equation (21a) becomes

[TABLE]

Likewise, we write $\theta^{k+1}=\theta^{k}+\Delta t\Omega^{k+1}$. By assuming $\Delta t\Omega\ll 1$, we obtain

[TABLE]

The tangent vector $\mathbf{r}^{k+1}_{s}$ at the filament base is opposite to the particle orientation $\mathbf{e}_{\mathrm{p}}$, namely $\mathbf{r}^{k+1}_{s}=\sin{\theta^{k+1}}\mathbf{e}_{y}-\cos{\theta^{k+1}}\mathbf{e}_{z}$. Using equation (63), the discretised form of this BC is,

[TABLE]

Combining equation (14a) and (23a), the discretised form for the force-free condition reads

[TABLE]

Combining equation (14b) and (23b) similarly, the torque-free condition reads

[TABLE]

We approximate $\theta^{k+1}$ by $\check{\theta}=2\theta^{k}-\theta^{k-1}$ and substitute it into equation (66), deriving the discretised torque-free condition

[TABLE]

Using the backward Euler scheme for $\mathcal{P}_{Q}$ and $\mathcal{P}_{3}$, and combining equations (63), we find the discretised governing equations for $\mathcal{P}_{Q}$ and $\mathcal{P}_{3}$,

[TABLE]

Integrating equations (18), (65), (66) and (68), with the four BCs equations (20a), (62) and (64) for $\mathbf{r}^{k+1}$ generates a linear system of size $3N+5$. Its solution corresponds to $[\mathbf{r}^{k+1}(s),U^{k+1}_{y},U^{k+1}_{z},\Omega^{k+1},\mathcal{P}^{k+1}_{Q},\mathcal{P}^{k+1}_{3}]$.

Appendix B Quincke rotation of a dielectric sphere jointed with a rigid

rod

Following Jones (1984), we derive the critical electric field required to trigger the electrohydrodynamic instability of a dielectric spherical particle grafted by a rigid rod. Let us first briefly reproduce the derivation of Jones (1984) for an individual particle and then extend it to our composite particle-rod configuration.

The electric torque exerted on a spherical particle of radius $A$ about its centre P is

[TABLE]

When the particle rotates about its centre at velocity $\Omega\mathbf{e}_{x}$, the hydrodynamic torque exerted on it is

[TABLE]

By using the torque-free condition $\boldsymbol{\Gamma}^{\mathrm{elec}}+\boldsymbol{\Gamma}^{\mathrm{hydro}}_{\mathrm{par}}=\mathbf{0}$, we derive

[TABLE]

Because the left-hand side of equation (71) is non-negative for a real value of $\Omega$, this condition gives us the criterion of the electrical field $E$ above which QR instability occurs,

[TABLE]

The rotational speed of the QR particle is known based on equation (71), so that its dimensionless value is

[TABLE]

where $\kappa=(R+2)/(S+2)$ as defined in equation (3).

Now we adapt the above derivation to the composite particle-rod system steadily rotating at velocity $\boldsymbol{\Omega}=\Omega\mathbf{e}_{x}$ about an off-centre pivot point V on the $y$-axis (see figure 17). We choose the particle centre P as the origin of the Cartesian coordinates. Using the local SBT, the force per unit length $-\mathbf{f}_{\mathrm{rod}}$ exerted by the fluid onto the rod at arclength $s$ is

[TABLE]

where $\xi_{\perp}=-8\pi\mu/c=-8\pi\mu/(1+2\log{\epsilon_{\mathrm{sl}}})$, and the total hydrodynamic force $\mathbf{F}_{\mathrm{rod}}$ exerted on the rod is

[TABLE]

Simultaneously, the hydrodynamic force exerted on the translating spherical particle is $\mathbf{F}_{\mathrm{par}}=6\pi\mu Ab\Omega\mathbf{e}_{z}$. Using the force-free condition $\mathbf{F}_{\mathrm{rod}}+\mathbf{F}_{\mathrm{par}}=\mathbf{0}$, we obtain

[TABLE]

The hydrodynamic torque $\boldsymbol{\Gamma}^{\mathrm{hydro}}_{\mathrm{rod}}$ exerted on the rod with respect to the particle centre P is

[TABLE]

Using the torque-free condition on the particle-rod system, $\boldsymbol{\Gamma}^{\mathrm{elec}}+\boldsymbol{\Gamma}^{\mathrm{hydro}}_{\mathrm{par}}+\boldsymbol{\Gamma}^{\mathrm{hydro}}_{\mathrm{rod}}=\mathbf{0}$, we obtain

[TABLE]

where

[TABLE]

Hence, the critical electrical field corresponding to the instability inception is

[TABLE]

The typical values of $\bar{\mathcal{E}}^{\mathrm{cri}}$ as a function of size ratio $\alpha=A/L$ for $\epsilon_{\mathrm{sl}}=0.01$ are provided in table 2.

Appendix C Two-dimensional polynomial equations for $\hat{\sigma}_{r}$ and

$W$

We substitute equations (47) and (48) into equation (46). We define $E^{\prime}=E^{2}-1$ and $\upsilon=\left(\frac{1}{-1-2\log{\epsilon_{\mathrm{sl}}}}\right)^{1/4}$ that leads to $\log\xi=z_{0}\upsilon W$. Consequently, the system of two-dimensional polynomial equations for $\hat{\sigma}_{r}$ and $W$ reads,

[TABLE]

and

[TABLE]

Appendix D LSA for the minimal model

We perform LSA for the minimal model. Similar to § 4.1, the bars over the dimensionless unknown variables are dropped, unless otherwise specified. The state variables $[\theta,\theta_{1},U_{y},U_{z},\mathcal{P}_{Q},\mathcal{P}_{3}]$ are decomposed into a base state $[\hat{\theta},\hat{\theta}_{1},\hat{U}_{y},\hat{U}_{z},\hat{\mathcal{P}}_{Q},\hat{\mathcal{P}}_{3}]$ and a perturbation state $[\theta^{\prime},\theta_{1}^{\prime},U_{y}^{\prime},U_{z}^{\prime},\mathcal{P}^{\prime}_{Q},\mathcal{P}^{\prime}_{3}]$, where $\hat{\theta}_{1}=\hat{U}_{y}=\hat{U}_{z}=0$, $\hat{\theta}$ is an arbitrary value, and $\hat{\mathcal{P}}_{Q}$ and $\hat{\mathcal{P}}_{3}$ are given by equation (26). By linearising equations (57) and (58) with respect to the base state, we derive the linear evolution equations for the perturbative state variables,

[TABLE]

Employing the normal-mode approach as in § 4, we substitute $\left[\theta^{\prime},\theta_{1}^{\prime},U_{y}^{\prime},U_{z}^{\prime},\mathcal{P}^{\prime}_{Q},\mathcal{P}^{\prime}_{3}\right]=\left[\Theta,\Theta_{1},\mathcal{Y},\mathcal{Z},\Phi,\Pi\right]\exp{(\sigma t)}$ into equation (D) and derive

[TABLE]

Setting the determinant of the operator matrix of equation (D) to zero, we find the non-zero solutions of the complex growth rate $\sigma$

[TABLE]

where

[TABLE]

We plot $\sigma_{r}$ and $\sigma_{i}$ in comparison with their numerical counterparts in figure 16.

Acknowledgments

We thank Drs. E. Han, L. Li, Y. Man and F. Yang, and Professors F. Gallaire, E. Nazockdast, O. S. Pak, B. Rallabandi and Y. N. Young for useful discussions. Prof. T. Götz is acknowledged for sharing with us his PhD thesis. We thank the anonymous referees for their insightful comments. L.Z. thanks the Swedish Research Council for a VR International Postdoc Grant (2015-06334). We thank the NSF for support via the Princeton University Material Research Science and Engineering Center (DMR-1420541). The computer time was provided by SNIC (Swedish National Infrastructure for Computing).