Stability of the interface of an isotropic active fluid

TL;DR

This paper analyzes the linear stability of isotropic active fluids in various geometries, revealing how activity influences interface stability and mode propagation, with implications for understanding active matter behaviors.

Contribution

It provides a comprehensive stability analysis of isotropic active fluids across different geometries using hydrodynamic theory, highlighting activity-induced instabilities and mode propagation.

Findings

Passive cases are stable in film and spherical droplet geometries.

Large activity induces instability in film and spherical droplet geometries.

Activity causes propagating damped or growing modes in all geometries.

Abstract

We study the linear stability of an isotropic active fluid in three different geometries: a film of active fluid on a rigid substrate, a cylindrical thread of fluid, and a spherical fluid droplet. The active fluid is modeled by the hydrodynamic theory of an active nematic liquid crystal in the isotropic phase. In each geometry, we calculate the growth rate of sinusoidal modes of deformation of the interface. There are two distinct branches of growth rates; at long wavelength, one corresponds to the deformation of the interface, and one corresponds to the evolution of the liquid crystalline degrees of freedom. The passive cases of the film and the spherical droplet are always stable. For these geometries, a sufficiently large activity leads to instability. Activity also leads to propagating damped or growing modes. The passive cylindrical thread is unstable for perturbations with wavelength longer than the circumference. A sufficiently large activity can make any wavelength unstable, and again leads to propagating damped or growing modes.

Equations

$${\mathcal{F}}_n=\int d^{3}x\left[\frac{A}{2}{Q}_{\alpha \beta }{Q}_{\alpha \beta }+\frac{B}{3}{Q}_{\alpha \beta }{Q}_{\beta \gamma }{Q}_{\gamma \alpha }+\frac{C}{4}(Q_{\alpha \beta }{Q}_{\alpha \beta }{)}^2\right],$$

$${\mathcal{F}}_i=\int \gamma \mathrm{d}S.$$

$$R_{\alpha \beta }={\partial }_t{Q}_{\alpha \beta }+\mathbf{v}\cdot \nabla Q_{\alpha \beta }+{\Omega }_{\alpha \gamma }{Q}_{\gamma \beta }-Q_{\alpha \gamma }{\Omega }_{\gamma \beta }.$$

$${\sigma }_{\alpha \beta }^v=2\eta E_{\alpha \beta }+2(\mu +{\mu }_1)R_{\alpha \beta }+a^{\prime }{Q}_{\alpha \beta },$$

$${\Phi }_{\alpha \beta }=2\mu E_{\alpha \beta }+\nu R_{\alpha \beta },$$

$$-A{Q}_{\alpha \beta }=2\mu E_{\alpha \beta }+\nu R_{\alpha \beta }.$$

$$Q_{\alpha \beta }=-\frac{2\mu }{A-\mathrm{i}\omega \nu }{E}_{\alpha \beta }.$$

$${\sigma }_{\alpha \beta }^v=2{\eta }_{\mathrm{eff}}{E}_{\alpha \beta },$$

$${\eta }_{\mathrm{eff}}$$

$$=$$

$${\partial }_{\beta }{\sigma }_{\alpha \beta }=0,$$

$${\eta }_{\mathrm{eff}}{\nabla }^2\text{v}-\nabla p=0.$$

$${\nabla }^4\mathbf{\psi }=0,$$

$${\mathbf{f}}_m=2\gamma H\mathbf{n},$$

$$({\sigma }_{nn}^{+}-{\sigma }_{nn}^{-})+2\gamma H$$

$$({\sigma }_{n\alpha }^{+}-{\sigma }_{n\alpha }^{-})+{\partial }_j{X}^{\alpha }$$

$$v_y(y=0)={\partial }_{t}h,$$

$${\sigma }_{xy}(y=0)=0,$$

$$\psi$$

$$F=\frac{\sinh 2kd-2kd}{\cosh 2kd+2k^2{d}^2+1}.$$

$$-{\sigma }_{yy}{|}_{y=0}+\gamma {\partial }_x^{2}h=0.$$

$$-\mathrm{i}\omega =-\frac{\gamma k}{2{\eta }_{\mathrm{eff}}(\omega )}F(kd).$$

$$-\mathrm{i}\omega$$

$$-\mathrm{i}\omega {\tau }_{\mathrm{lc}}^{\prime }=\frac{a-1}{2}-\frac{kdF{\tau }_{\mathrm{lc}}}{4{\tau }_{\mathrm{s}}}\pm \sqrt{{\left(\frac{a-1}{2}-\frac{kdF{\tau }_{\mathrm{lc}}}{4{\tau }_{\mathrm{s}}}\right)}^2-\frac{kdF{\tau }_{\mathrm{lc}}^{\prime }}{2{\tau }_{\mathrm{s}}}},$$

$$a_{\mathrm{c}}(k)=1+\frac{1}{2}kdF(kd){\tau }_{\mathrm{lc}}/{\tau }_{\mathrm{s}}.$$

$$a_{\pm }=1+\frac{1}{2}kdF(kd){\tau }_{\mathrm{lc}}/{\tau }_{\mathrm{s}}\pm \sqrt{\frac{1}{2}kdF(kd){\tau }_{\mathrm{lc}}^{\prime }/{\tau }_{\mathrm{s}}}.$$

$$D^2\Psi =0,$$

$$\frac{k}{R}\Psi (\rho =R)=-\omega ϵ.$$

$$\psi =ϵ\omega {\mathrm{e}}^{\mathrm{i}kx}\left[\frac{{\rho }^2{I}_0(k\rho )}{I_1(kR)}-\frac{kR{I}_0(kR)+I_1(kR)}{k{I}_1^2(kR)}\rho I_1(k\rho )\right].$$

$$-{\sigma }_{\rho \rho }{|}_{\rho =R}+2\gamma H=0.$$

$$H=-\frac{1}{2}\left[\frac{1}{R}+ϵ\left(k^2-\frac{1}{R^2}\right){\mathrm{e}}^{\mathrm{i}kx}\right],$$

Full text

Stability of the interface of an isotropic active fluid

**Harsh Sonia, Wan Luoa, Robert A. Pelcovitsb, and Thomas Powersa,b

**

We study the linear stability of an isotropic active fluid in three different geometries: a film of active fluid on a rigid substrate, a cylindrical thread of fluid, and a spherical fluid droplet. The active fluid is modeled by the hydrodynamic theory of an active nematic liquid crystal in the isotropic phase. In each geometry, we calculate the growth rate of sinusoidal modes of deformation of the interface. There are two distinct branches of growth rates; at long wavelength, one corresponds to the deformation of the interface, and one corresponds to the evolution of the liquid crystalline degrees of freedom. The passive cases of the film and the spherical droplet are always stable. For these geometries, a sufficiently large activity leads to instability. Activity also leads to propagating damped or growing modes. The passive cylindrical thread is unstable for perturbations with wavelength longer than the circumference. A sufficiently large activity can make any wavelength unstable, and again leads to propagating damped or growing modes.

††footnotetext: a* School of Engineering, Brown University, Providence, RI 02912, USA. ††footnotetext: b Department of Physics, Brown University, Providence, RI 02912, USA. *

1 Introduction

Active fluids are energized locally by motorized microscopic active particles such as kinesin-driven microtubules 1 and myosin-actin complexes 2. Therefore, their dynamics occur out of thermal equilibrium 3, 4. Hydrodynamic instabilities of both polar and nematic active fluids have been studied using hydrodynamic theories and simulations for bulk 5, 6, 7 as well as for confined fluids 8, 9, 10, 11, 12, 13, 14. In this paper, we consider the instabilities of active nematic fluids in the isotropic phase confined by an interface. The damping of a capillary wave on a flat interface between two passive viscous fluids is well-understood 15. Likewise, theoretical studies of interfacial instabilities like the Rayleigh-Plateau capillary instability and Rayleigh-Taylor interface instability have been carried out for passive fluids 16 including complex fluids such as polymer solutions 17, 18 and liquid crystals 19, 20. Less work has been done on interfacial instabilities in active fluids. It is natural to expect that the instabilities that occur in bulk active fluids can destabilize an otherwise stable interface, or make an already unstable interface more unstable. Work to date includes a study by Yang and Wang 21 of the Rayleigh-Plateau capillary instability of a thread of active polar fluid in the ordered state surrounded by a passive Newtonian fluid, a study by Whitfield and Hawkins 22 of the instability of an spherical droplet of active polar fluid in the ordered state, and an analysis by Gao and Li 23 of a self-driven droplet of an active nematic fluid. Also, Patteson et al. 24 studied the propagation of active-passive interfaces in bacterial swarms. Recently, Maitra et al. 25 explored the dynamics of an active membrane in an active polar medium, Mietke et al. studied the instabilities of an active membrane in a passive fluid 26, and V. Soni et al. studied the surface dynamics of an active colloidal chiral fluid 27. Here, we focus on linear stability analyses of active nematic fluids in the isotropic phase in with flat, cylindrical, or spherical interfaces. Our focus on the isotropic phase is motivated by recent experiments on active matter that show large regions in which the nematic order is small 28. Our work is distinct from the other theoretical work just mentioned on interfacial instabilities in active fluids because we consider the active nematic to initially be in the isotropic state instead of the ordered state.

We model the isotropic phase of an active nematic fluid by adding activity to de Gennes’ hydrodynamic model 29, 30, 31, 32 for the isotropic phase of a passive nematic. This model is appropriate for ‘shakers’ rather than ‘movers’ suspended in a liquid. The model shows that in the linear regime the isotropic active nematic fluid behaves like a viscoelastic fluid, with the viscosity and viscoelasticity growing large as the isotropic-nematic transition is approached 33. However, our results for the stability of interfaces are qualitatively different from the passive viscoelastic fluid case due to the orientational degrees of freedom. We work in the limit of low Reynolds number, where viscous effects dominate inertial effects. For a passive fluid, deformations of a surface or spherical surface always relax, whereas a cylindrical thread is unstable to peristaltic deformations of sufficiently long wavelength. When the fluid is active, deformations of the surfaces in all three cases can be unstable. The instability of a bulk active isotropic fluid drives the instability of the flat and spherical surface, and enhances the Rayleigh-Plateau capillary instability of a cylinder.

Our key results are as follows. In the all cases we consider, the coupled dynamics of the interface of the fluid and the nematic order parameter leads to two modes, with damped or growing propagating waves found for a sufficiently large dimensionless activity. Likewise, surface tension makes it harder for activity to destabilize an active film or an active fluid confined by a spherical interface, as compared to the unconfined case. For the cylindrical thread of radius $R$, harmonic perturbations of wavenumber $k<1/R$ are always unstable, just as in the passive case. Perturbations with $k>1/R$ become unstable above a critical activity increasing with $k$ and the surface tension of the interface.

The remainder of the paper is organized as follows. In section 2, we introduce a hydrodynamic model for an active nematic fluid in the isotropic state. In section 3, we use this model to study the linear stability of a film bound by an interface. Next , in section 4, we consider the stability of a thread of active fluid bound by either an interface. Finally, in section 5 analyze the stability of a spherical drop of active fluid. We offer concluding remarks in Section 6. Section 7, the Appendix, contains additional details relevant to Section 3.

2 Model

The total free energy of an active isotropic nematic fluid with an interface is $\mathcal{F}=\mathcal{F}_{n}+\mathcal{F}_{i}$, where $\mathcal{F}_{n}$ is the free energy of the nematic fluid, and $\mathcal{F}_{i}$ is the energy of the interface. Denoting the nematic order parameter field by $Q_{\alpha\beta}$, the nematic free energy is

[TABLE]

where we sum over repeated indices $\alpha,\beta,\dots$ which run over the three spatial coordinates. We consider the isotropic phase, for which $A>0$. In this case, Frank elasticity can be neglected as long as we are not too near the nematic transition.

The interface energy is given by

[TABLE]

where $\gamma$ is the interfacial tension and $\mathrm{d}S$ is the element of area.

We use de Gennes’ hydrodynamic model 30, 31, 32 of the isotropic phase of a passive nematic fluid of uniform concentration, suitably modified 29 to account for activity. In terms of the fluid velocity field $v_{\alpha}$, the strain rate and the vorticity tensors are given by $E_{\alpha\beta}=(\partial_{\alpha}v_{\beta}+\partial_{\beta}v_{\alpha})/2$ and $\Omega_{\alpha\beta}=(\partial_{\alpha}v_{\beta}-\partial_{\beta}v_{\alpha})/2$, respectively, where $\alpha,\beta=x,y,z$. The rate of change $R_{\alpha\beta}$ of the nematic order parameter $Q_{\alpha\beta}$ relative to the local background fluid is defined as

[TABLE]

Then, the viscous stress $\sigma^{v}_{\alpha\beta}$ and equation of motion for the nematic order parameter $Q_{\alpha\beta}$ are given by 29

[TABLE]

where $\eta$ and $\nu$ are the shear and rotational viscosities, respectively, $\mu$ couples shear and nematic alignment, $a^{\prime}$ and $\mu_{1}$ are activity parameters, and $\Phi_{\alpha\beta}$ is the molecular field defined as $\Phi_{\alpha\beta}\equiv-\delta\mathcal{F}/\delta Q_{\alpha\beta}$. From Eq. (5) we can see that for the case of small shear rate and steady state, the principal axes of the order parameter align with the principal axes of the strain, with the case $\mu<0$ corresponding to the way that prolate particles align in shear, and the case $\mu>0$ corresponding to the way oblate particles align in shear.

In passive fluids, $a^{\prime}=0$ and $\mu_{1}=0$. In that case, the Onsager reciprocal relations 34 are obeyed and the positive entropy production rate leads to the relation $\eta\nu-2\mu^{2}>0$. The active term $aQ_{\alpha\beta}$ appearing in Eq. (4) accounts for the stress due to the force dipoles associated with the active particles 5, 29 with $a^{\prime}>0$ for contractile particles and $a^{\prime}<0$ for extensile particles. Since time reversal symmetry and the Onsager relations are violated in active fluids when we do not keep track of the chemical reactions in the theory, the active term $\mu_{1}R_{\alpha\beta}$ is allowed in Eq. (4). Various other approaches to the active matter equations also violate the Onsager reciprocal relation 5, 10, 12. In our entire analysis, we assume that $\mu_{1}$ is sufficiently small such that $\eta\nu-2\mu(\mu+\mu_{1})>0$. Since we only study linear stability of the state with no order and no flow, we are justified in disregarding terms of higher order than quadratic in the order parameter. Thus, $\Phi_{\alpha\beta}\approx-AQ_{\alpha\beta}$ (with $A>0$ in the isotropic phase) and Eq. (5) takes the form

[TABLE]

Likewise, we ignore the higher order terms in $R_{\alpha\beta}$; thus, $R_{\alpha\beta}\approx\dot{Q}_{\alpha\beta}$. Our linearized equations are equivalent to the apolar case of the linearized equations of active matter that have appeared previously 33, 35, 36, 37 when we set $\mu_{1}=0$; also we absorb a possible active term proportional to $Q_{\alpha\beta}$ in $\Phi_{\alpha\beta}$ in Eq. (5).

Assuming that $v_{\alpha}$ and $Q_{\alpha\beta}$ are proportional to $\exp(-\mathrm{i}\omega t)$, where the real part of $-\mathrm{i}\omega$ is the growth rate of the perturbations, we find using Eq. (6) that

[TABLE]

Using Eq. (4), the viscous stress $\sigma^{v}_{\alpha\beta}$ is given by

[TABLE]

where

[TABLE]

Here, we have defined the dimensionless activity by $a=a^{\prime}\mu/\eta A$, and the relaxation times by $\tau_{\mathrm{lc}}=\nu/A$ and $\tau_{\mathrm{lc}}^{\prime}=\tau_{\mathrm{lc}}[1-2\mu(\mu+\mu_{1})/(\nu\eta)]$. Since we are assuming that $\eta\nu-2\mu(\mu+\mu_{1})>0$, the isotropic phase of an infinite active nematic fluid is unstable against shear flow and local ordering 29 when $a>1$. We will see in our instability analyses for the various geometries that the critical values of the dimensionless activity correspond to a negative effective shear viscosity, i.e. $a\geq 1$. Note that in the oblate particle case of $\mu>0$, the critical value of the activity $a^{\prime}$ is positive, meaning that significantly active contractile (puller) particles lead to instablity. For the prolate particle case of $\mu<0$, the critical activity is negative, meaning that sufficiently active extensile (pusher) particles lead to instability.

It is apparent from the above equation that the effective viscosity of the fluid $\eta_{\mathrm{eff}}$ depends on the growth rate $-\mathrm{i}\omega$; in other words, the fluid behaves like a viscoelastic fluid due to the presence of the nematic molecules. At the special value of dimensionless activity $a=2\mu(\mu+\mu_{1})/\nu\eta=1-\tau_{\mathrm{lc}}^{\prime}/\tau_{\mathrm{lc}}$, the effective viscosity $\eta_{\mathrm{eff}}$ is independent of $\omega$, and the fluid behaves like a Newtonian fluid with shear viscosity $\eta_{\mathrm{eff}}=\eta-2\mu(\mu+\mu_{1})/\nu$. We will see below that at this special value of the activity the growth rate is that of a passive fluid.

Since we ignore the inertia of the fluid, the force balance equation is given by

[TABLE]

with $\sigma_{\mathrm{\alpha\beta}}=-p\delta_{\alpha\beta}+\sigma^{v}_{\alpha\beta}$. The pressure $p$ is the pressure arising from the incompressibility condition $\bm{\nabla}\cdot\textbf{v}=0$. We have disregarded the Ericksen stress $\sigma^{\mathrm{e}}_{\alpha\beta}=\mathcal{F}\delta_{\alpha\beta}-\partial{\mathcal{F}}/\partial(\partial_{\beta}Q_{\mu\nu})\partial_{\alpha}Q_{\mu\nu}$ 30, 38 since it is at least quadratic order in $Q_{\alpha\beta}$. Then, using the incompressibility condition $\bm{\nabla}\cdot\textbf{v}=0$, the linearized Eq. (11) can be simplified to

[TABLE]

The incompressibility condition is imposed by representing v as the curl of a stream function $\bm{\psi}$ i.e. $\textbf{v}=\bm{\nabla}\times\bm{\psi}$. For simplicity, we choose the form of $\bm{\psi}$ such that $\bm{\nabla}\cdot\bm{\psi}=0$. Taking the curl of Eq. (12) yields

[TABLE]

where $\nabla^{4}$ is the square of the Laplacian operator in three dimensions. We solve the above equation with the boundary conditions appropriate to the geometry at hand and calculate the forces on the interface due to the fluid.

To describe the force per unit area acting on at the surface, we need to parametrize the surface as $\mathbf{X}(u^{1},u^{2})$, with coordinates $u^{1}$ and $u^{2}$. Due to the free energy associated with the surface [see Eq. (2)], the force per unit area acting on the surface is given by 39

[TABLE]

where $\mathbf{n}$ is the outward normal. Note that our convention is that $H$ is negative for a sphere or a cylinder. Since we disregard the inertia of the surface, the force balance equation at the surface reads

[TABLE]

where $\sigma^{\pm}_{nn}=n^{\alpha}\sigma^{\pm}_{\alpha\beta}n^{\beta}$ and $\sigma_{n\alpha}=n^{\beta}\sigma_{\beta\alpha}$, with the plus and minus denoting the stress exerted on the interface from the $\mathbf{n}$ and $-\mathbf{n}$ sides, respectively.

We close this section with estimates of the magnitudes of the liquid crystal relaxation time $\tau_{\mathrm{lc}}$ and the characteristic time scales for a film with interfacial tension or bending stiffness. A crude dimensional analysis estimate for $\tau_{\mathrm{lc}}=\nu/A$ is to suppose $\nu\approx\eta$, and to take $A=k_{\mathrm{B}}T/\ell^{3}$, where $k_{\mathrm{B}}T$ is thermal energy and $\ell$ is the length of the active particles. Using the viscosity of water, $\eta\approx 10^{-3}\,$N-s/m2, and $\ell\approx 10\,\mu$m leads to $\tau_{\mathrm{lc}}\approx 300\,$s. If the rods are $1\,\mu$m in length, then $\tau_{\mathrm{lc}}\approx 0.3\,$s. However, since we are considering an active system, it is reasonable to suppose that $A$ is not determined by thermal energy, and that $A$, and the liquid crystal relaxation rate may be much bigger. For a film of thickness $d\approx 1\,$mm and for the air-water surface tension $\gamma\approx 70\times 10^{-3}\,$N/m, the characteristic surface-tension driven relaxation time is $\tau_{\mathrm{s}}=\eta d/\gamma\approx 0.1\,$ms. Thus, we expect the film relaxation time to be much shorter than the liquid crystal relaxation time, and we will focus on this limit. However, due to our uncertainty about the value of $A$, and also to show some of the range of possible phenomena, we also consider the case of $\tau_{\mathrm{lc}}\approx\tau_{\mathrm{s}}$.

3 Instability of an active fluid film

In this section, we study the instability of a flat interface of an active nematic fluid in its isotropic phase. The fluid is a film of thickness $d$ atop a solid substrate, with air above the film (Fig. 1).

We consider an air-fluid interface with constant uniform surface tension $\gamma$, and no bending stiffness. A film of passive fluid is always stable to sinusoidal perturbation, since the perturbation increases the surface area. Thus, the instability we study in this section arises from the activity of the fluid.

The surface, which lies in the $zx$ plane in its unperturbed state, is subject to a transverse perturbation which is the real part of $h=\epsilon(t)\exp(\mathrm{i}kx)$, as shown in Fig. 1. We assume that $\epsilon\propto\exp(-\mathrm{i}\omega t)$.

The stream function is given by $\bm{\psi}=\psi\hat{\mathbf{z}}$, with $\psi$ a biharmonic function. For small deflections $kh\ll 1$, the kinematic condition takes the form

[TABLE]

where $v_{y}=-\partial_{x}\psi$. This condition, along with the conditions of zero tangential stress at the interface,

[TABLE]

and vanishing flow at $y=-d$, leads to

[TABLE]

where

[TABLE]

For small deflections the mean curvature is $H\approx-\partial^{2}_{x}h/2$, and the force balance equation on the interface becomes

[TABLE]

The stress component $\sigma_{yy}$ can be found by calculating the pressure from the $x$-component of the force balance equation (11). Once $\sigma_{yy}$ is calculated, we use normal stress balance (21) at $y=0$ to obtain the characteristic equation

[TABLE]

In the passive Newtonian case with $a=0$ and with no coupling between the fluid and the liquid crystalline degrees of freedom, i.e. $\mu=\mu_{1}=0$, the growth rate has two branches that cross, one corresponding to the negative growth rate of a Newtonian film 40, with characteristic time scale $\tau_{\mathrm{s}}=\eta d/\gamma$,

[TABLE]

and one corresponding to the liquid crystalline relaxation rate, $-\mathrm{i}\omega=-1/\tau_{\mathrm{lc}}=-A/\nu$. When $\mu$ (or $\mu_{1}$) is nonzero and $a=0$, the growth rate curves repel each other instead of crossing, as in Fig. 2, upper left panel.

The active case is like the case of a passive viscoelastic fluid 40, for which the effective shear viscosity depends on $\omega$, and we must solve Eq. (22) for $\omega$ as a function of $k$, which yields

[TABLE]

where $F$ is given by Eq. (20). As the activity increases, the splitting between the two growth rate curves decreases, until the value of $a=2\mu\left(\mu+\mu_{1}\right)/\eta\nu$ is reached. At this special value of activity, $\eta_{\mathrm{eff}}$ is independent of $\omega$, and the branches of the growth rates cross as they do in the case of $a=0$ and $\mu=0$ (Fig. 2, upper right panel). As the activity increases further, the real branches collapse into one branch for a range of wavevector, and the imaginary parts of the growth rate become nonzero in this same range (Fig. 2, lower left panel). The critical activity $a=1$ corresponds to the point at which the effective shear viscosity vanishes. When $a>1$, one of the branches of the real part of the growth rate becomes positive, and the system is unstable for sufficiently long wavelengths (Fig. 2, lower right panel). The critical activity $a_{\mathrm{c}}(k)$ at which the mode $k$ is marginally stable is found by demanding that $\mathrm{Re}(-\mathrm{i}\omega)=0$:

[TABLE]

Since $\tau_{\mathrm{lc}}/\tau_{\mathrm{s}}=(\nu\gamma)/(\eta Ad)$, interfacial tension tends to suppress the instability for nonzero $k$. But even if $\tau_{\mathrm{s}}\ll\tau_{\mathrm{lc}}$, the longest wavelengths are always unstable. In this limit, the two branches of the uncoupled passive case cross when $kd\sim(\tau_{\mathrm{s}}/\tau_{\mathrm{lc}})^{1/4}$, which is why we plot the growth rates vs. $kd(\tau_{\mathrm{lc}}/\tau_{\mathrm{s}})^{1/4}$ in Fig. 2. The shapes of the real and imaginary parts of the growth rate curves for $\tau_{\mathrm{lc}}\approx\tau_{\mathrm{s}}$ and for $\tau_{\mathrm{lc}}\ll\tau_{\mathrm{s}}$ are qualitatively similar to the case of $\tau_{\mathrm{s}}\ll\tau_{\mathrm{lc}}$, with the main difference being that the band of unstable modes reaches further into the regime of short wavelength as $\tau_{\mathrm{s}}$ increases relative to $\tau_{\mathrm{lc}}$ (See Figs. 10 and 11 in the appendix).

Note that the growth rate $-\mathrm{i}\omega$ always has an imaginary part when $a$ is sufficiently near $a_{\mathrm{c}}(k)$; when a mode is unstable with a sufficiently small growth rate, it also propagates. Propagating modes are found when $a_{-}\leq a\leq a_{+}$, where

[TABLE]

Also, there are no propagating modes without the interface, since $a_{+}-a_{-}\propto\sqrt{\tau^{\prime}_{\mathrm{lc}}/\tau_{\mathrm{s}}}\propto\sqrt{\gamma}$.

Figure 3 shows when the interface is stable as a function of scaled dimensionless wavenumber and dimensionless activity for the case of $\tau_{\mathrm{lc}}/\tau_{\mathrm{s}}\gg 1$. The system is always stable for $a<1$; as $a$ is increased beyond $a=1$, an increasingly large band of very long wavelength modes are unstable. Growing and decaying modes with a sufficiently small growth rate [between the dashed lines in Fig. 3, which are given by Eq. (27)] are also propagating. As $\tau_{\mathrm{lc}}/\tau_{\mathrm{s}}$ decreases, the band of unstable modes is limited to shorter and shorter wavelengths. To sum up, the time scale $\tau_{\mathrm{lc}}$ controls the rate of growth or decay of the modes, and the time scale $\tau_{\mathrm{s}}$ determines which modes become unstable. Since $a_{\pm}$ depends on $\tau_{\mathrm{s}}$, the velocity of propagation $\mathrm{Re}(\omega)/k$ is determined by $d/\tau_{\mathrm{s}}$.

4 Rayleigh-Plateau capillary instability

A fluid thread breaks into drops because perturbations of sufficiently long wavelength lower the area of the surface, and thus the energy. This instability is known as the Rayleigh-Plateau capillary instability 41, 42. In this section, we study how the presence of active nematic molecules in the liquid affects the Rayleigh-Plateau capillary instability. For simplicity, we disregard the outer fluid. While this approximation was natural in our study of the stability of a flat interface between air and an active fluid, it seems less natural for a thread of active fluid, since the thread must be supported by some surrounding fluid if it is not a jet. However, unlike the passive case of a stationary cylindrical interface 16, accounting for the viscosity contrast leads to a complicated characteristic equation for the growth rate of the interface of an active thread. To avoid this complication and illustrate the essential physics, we assume the outer fluid is of sufficiently small viscosity that we may disregard it.

We consider a cylindrical fluid thread of initial radius $R$, subject to an axisymmetric harmonic perturbation of wavenumber $k$ along the $x$ direction (see Fig. 4). The cylindrical coordinates are $(\rho,\theta,x)$. Initially the fluid is at rest, with a uniform pressure $p=\gamma/R$. The radius of the perturbed thread is given by the real part of $h(x,t)=R+\epsilon(t)\exp(\mathrm{i}kx)$, with $\epsilon k\ll 1$. For an axisymmetric flow, we follow Happel and Brenner 43 and define the stream function via $\bm{\psi}=-(\psi/\rho)\hat{\bm{\theta}}$. The stream function $\psi$ is related to velocity by $v_{\rho}=(1/\rho)\partial_{x}\psi$ and $v_{x}=-(1/\rho)/\partial_{\rho}\psi$. If we choose $\bm{\psi}=\Psi(\rho)\exp(\mathrm{i}kx)\hat{\bm{\theta}}$, then Eq. (13) in cylindrical coordinates reduces to

[TABLE]

where 16 $D\equiv\partial_{\rho}^{2}-(1/\rho)\partial_{\rho}-k^{2}$. The linearized kinematic condition at the interface, $\bm{\nabla}\times\bm{\psi}=\partial_{t}{h}\hat{\bm{\rho}}$, leads to

[TABLE]

The kinematic boundary condition and the condition of zero tangential stress, $\sigma_{x\rho}|_{\rho=R}=0$, along with the condition of regularity at $\rho=0$, leads to the solution

[TABLE]

where $I_{0}$ and $I_{1}$ are the Bessel functions of first kind. The growth rate is determined by the normal force balance equation,

[TABLE]

The pressure may be found from the $x$-component of the Stokes equation, Eq. (12); with this pressure and the velocity field we may calculate $\sigma_{\rho\rho}=-p+2\eta_{\mathrm{eff}}\partial_{\rho}v_{\rho}$ and use the mean curvature expanded 44 to linear order in $\epsilon$,

[TABLE]

in Eq. (31) to find

[TABLE]

where

[TABLE]

When $\eta_{\mathrm{eff}}(\omega)=\eta$, the growth rate of Eq. (33) is precisely that of a thread of a passive Newtonian viscous fluid thread 45. Since the characteristic equation (33) for the cylinder is of a similar form as the characteristic equation (22) for the planar surface, the growth rate is given by Eq. (25) with $F$ replaced by $-G/k$ and $d$ replaced by $R$. (Note that in this section $\tau_{\mathrm{s}}=\eta/(\gamma R)$.) Figure 5 shows the growth rate vs. dimensionless wavevector $kR$ for the case of $\tau_{\mathrm{lc}}\gg\tau_{\mathrm{s}}$. In this case, the growth rate is almost exactly the same as the classical result for a passive Newtonian fluid. The only dependence on activity or liquid crystalline parameters arises in the region near $kR=1$ where the real part of the growth rate vanishes. This fact can be seen by expanding the growth rate for small $\tau_{\mathrm{s}}/\tau_{\mathrm{lc}}$; away from the region where $G\ll 1$, we have

[TABLE]

The effects of activity become apparent when the liquid crystal relaxation time is comparable to the film relaxation time, $\tau_{\mathrm{lc}}\sim\tau_{\mathrm{s}}$. The growth rate for several different dimensionless activities is shown in Fig. 6. In this case, the behavior of the growth rate with respect to activity is similar to behavior of the growth rate for a flat interface (compare with Fig. 2). The passive cylindrical thread is always unstable for modes with $kR<1$. Likewise, in the active case, modes with $kR<1$ are always unstable. Once $a>1$, modes with a wavenumber greater than $1/R$ can also be unstable; in particular, $\mathrm{Re}(-\mathrm{i}\omega)=0$ when

[TABLE]

Propagating modes are found when $a_{-}<a<a_{+}$, where

[TABLE]

Note that propagation only occurs when $kR>1$, i.e. $G(k)<0$. Figure 7 is the stablity diagram for the case of $\tau_{\mathrm{lc}}=\tau_{\mathrm{s}}$.

5 Instability of a spherical active droplet

A cylinder of active fluid is unstable, and breaks up into spherical droplets. A spherical droplet of a Newtonian fluid is always stable against surface tension since the spherical shape minimizes the surface energy. However, a spherical droplet of active fluid might go unstable due to activity. Here we carry out a linear stability analysis for a droplet of active nematic fluid in the isotropic phase (see Fig. 8). We assume that the spherical droplet of radius $R$ is subject to spherical harmonic perturbations such that the surface of the perturbed drop can be represented by $\mathbf{X}(\theta,\phi)=(R+\epsilon(t)Y^{m}_{l}(\theta,\phi))\hat{\bm{r}}$ with $\epsilon\ll R$. We choose the following form of the stream function $\bm{\psi}$ to enforce the condition $\bm{\nabla}\cdot\bm{\psi}=0$:

[TABLE]

where $Y^{m}_{l}(\theta,\phi)$ is a spherical harmonic. Inserting this stream function in the Stokes equations, we find that the function $v(r)$ obeys

[TABLE]

where

[TABLE]

The boundary conditions on the interface are the linearized kinematic condition, $\bm{\nabla}\times\bm{\psi}=\partial_{t}\textbf{X}$, and the linearized zero shear stress condition:

[TABLE]

The solution of Eq. (40) with the above boundary conditions is given by

[TABLE]

With this solution, we get the following expression for $\sigma_{rr}$ after integrating the $r$-component of Eq. (13) with respect to $r$:

[TABLE]

where

[TABLE]

In the unperturbed state, the surface tension leads to a constant pressure $C$ via the Young-Laplace law. Since we suppose that there is no fluid outside the drop, the force balance equation at the surface of the drop (in the limit $\epsilon\ll R$) is given by (see Eq. (16))

[TABLE]

The mean curvature $H$ is given to first order in $\epsilon$ by 44,

[TABLE]

We see from Eqs. (46) and (49) that, for the $l=1$ mode, there are no changes in $\sigma_{rr}$ or the Laplace pressure $2\gamma H$ due to the perturbation, because to leading order, the $l=1$ mode is equivalent to the displacement of the droplet along the $z$ direction (see Fig. 8). Therefore, we consider modes with $l>1$. From Eq. (46), (48) and (49), we find that $C=-2\gamma/R$ and

[TABLE]

When $\eta$ is independent of $\omega$, this result is precisely the relaxation rate for perturbations of a sphere with surface tension in the limit that viscosity dominates inertia 46, 47.

Equation (50) is quadratic in $\omega$, and the real parts of its two roots represent growth rates of the perturbation. The critical dimensionless activity $a_{c}(l)$ for the $l$th harmonic perturbation calculated is given by

[TABLE]

Since the smallest value of $l$ is 2, the critical value of the dimensionless activity above which droplet becomes unstable is given by

[TABLE]

Therefore, critical dimensionless activity for a spherical droplet $a_{c}(l=2)$ is larger than its value for the unconfined fluid. Also, $a_{c}(l=2)$ decreases with $R$: smaller active droplets are more stable. Fig. 9 shows that $a_{c}(l)$ increases almost linearly with $l$.

6 Discussion and Conclusion

In this paper we have studied the effect of activity on the stability of flat, cylindrical, and spherical interfaces. In all cases, the bulk instability of the active fluid, which is characterized by a vanishing effective shear viscosity, leads to spontaneous shear flows that can destabilize an interface that would be stable in the case of a passive fluid. Furthermore, all three geometries showed oscillatory behavior at suitably large activity, corresponding to propagating damped or growing modes. The presence of propagating modes (damped or growing) at zero Reynolds number is qualitatively different from the passive fluid case, where no propagation is seen at zero Reynolds number. The propagating modes in our linear stability analysis may be the seed for propagating modes at large amplitude, as seen in numerical calculations of active membranes 26. We made several approximations in this paper to make our calculation tractable. We neglected the Frank elasticity, which meant that the base state that we expanded about is uniform, $Q_{\alpha\beta}=0$. If we had included Frank elasticity, we would have to specify anchoring conditions for $Q_{\alpha\beta}$. For the case of planar or homeotropic anchoring, the base state would be nonuniform, and its stability would be more difficult to analyze by the technique we employ. The case of a zero-torque anchoring condition would lead to a uniform base state, but it would still make our calculation more complicated since we would not be able to eliminate $Q_{\alpha\beta}$ by simply solving an algebraic equation, and we would not be able to lump all the liquid-crystalline and active effects into the effective frequency-dependent viscosity $\eta_{\mathrm{eff}}(\omega)$. It would be interesting to generalize our calculations to include Frank elasticity, since it has been shown that Frank elasticity (or equivalently rotational diffusion in the work of Woodhouse and Goldstein) leads to spontaneous flow even for undeformed confining surfaces 10. A second major simplification is our neglect of the outer fluid. Because we neglected the viscosity of the outer fluid, we only had to solve a quadratic equation to find the branches of the growth rate. Including the outer fluid is more realistic, and it will lead to a more complicated characteristic equation, and more branches. Also, if we use the thermal energy scale to estimates the material parameters (questionable in a active system), we are led to $\tau_{\mathrm{lc}}\gg\tau_{\mathrm{s}}$, which makes the interesting activity-driven phenomena such as instability and oscillation occur at long wavelength in the case of the flat film, but only in a narrow regime near $kR\approx 1$ in the case of the cylindrical thread. When the viscosity of the outer fluid is accounted for, the growth rate of the passive cylindrical thread vanishes 16 at $k=0$, which will also lead to interesting activity-driven behavior at long wavelength in the cylinder. Finally, all of the calculations we did for interfaces could be modified to apply to the case of an active fluid bound by a membrane, which could be more relevant for biological phenomena.

7 Appendix

In this appendix we display more plots of the growth rate and the stability diagram for the case of the film of thickness $d$ (Section 3). Fig. 10 shows the real and imaginary parts of the growth rate for $\tau_{\mathrm{lc}}=\tau_{\mathrm{s}}$, whereas Fig. 11 shows the same quantities for the case of $\tau_{\mathrm{s}}/\tau_{\mathrm{lc}}\gg 1$. In all case, the shape of the curves is qualitatively similar, but the scale of wavevectors where the instability and oscillations changes, with the instability and oscillations occurring when $kd\sim(\tau_{\mathrm{s}}/\tau_{\mathrm{lc}})^{1/4}$ when $\tau_{\mathrm{s}}/\tau_{\mathrm{lc}}\ll 1$, when $kd\sim 1$ when $\tau_{\mathrm{s}}/\tau_{\mathrm{lc}}\sim 1$, and when $kd\sim\tau_{\mathrm{s}}/\tau_{\mathrm{lc}}$ when $\tau_{\mathrm{s}}/\tau_{\mathrm{lc}}\gg 1$. Figure 12 shows the stability diagram for $\tau_{\mathrm{s}}=\tau_{\mathrm{lc}}$ (upper panel) and $\tau_{\mathrm{s}}\gg\tau_{\mathrm{lc}}$ (lower panel).

8 Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported in part by National Science Foundation Grant Nos. CBET-1437195 (TRP) and National Science Foundation Grant No. MRSEC-1420382 (RAP and TRP). We are grateful to Dan Blair, Kenny Breuer, and Ian Wong for helpful discussions.