Instability of a thin film of chemotactic active suspension

Abstract

In this paper, we analyse the instability of a thin-film suspension of micro-swimmers subjected to an attractant gradient. The imposed attractant gradient introduces a preferential swimming direction for the swimmers, resulting in an anisotropic orientation field. By modelling the swimmers as force dipoles, the study by Kasyap & Koch (Phys. Rev. Lett., vol. 108, issue 3, 2012, 038101, J. Fluid Mech., vol. 741, 2014, pp. 619–657) found an instability arising from the coupling between the active stress associated with an anisotropic orientation field and perturbations to the swimmer density field. In this work we begin by presenting a stability analysis that calculates the modification of this instability in the presence of an interface. We then detail the presence of a new mode of instability that arises solely from the deformation of the interface mediated by the active stress in the suspension. The resulting hydrodynamic instabilities in the system are observed to be sensitive to the direction of the attractant gradient relative to the interface. Furthermore, we show that the coupling between the two modes involving a jump in interfacial viscous stresses and normal stress differences within the suspension allows for an instability to manifest in a suspension of chemotactic pullers, previously thought to be unconditionally stable. We then use a long-wave theory to write down a set of nonlinear equations governing the film evolution. The numerical solution of the system using the long-wave theory aids in explaining the mechanism associated with the instability in addition to validating the predictions from the linear stability theory.

1. Introduction

Theoretical and experimental investigations have revealed the existence of flow instabilities in active suspensions of fluids (Ramaswamy 2010; Koch & Subramanian 2011). Unlike passive particles in fluids, active suspensions contain motile organisms capable of individual motion. The cell body of each swimmer contains one or more flagella or cilia which act like motors and churn the surrounding fluid to generate a propulsive force. The Reynolds number for a swimmer such as Escherichia coli, which is approximately $7\ \mathrm {\mu } {\rm m}$ long, in an aqueous medium with a swimming speed of $30 \ \mathrm {\mu } {\rm m}\ {\rm s}^{-1}$ is ${\sim }{O}(10^{-4})$. The negligible magnitude of inertial forces associated with the swimmer requires that the propulsive force generated by the rotating flagella or cilia be balanced by the drag on the cell body (Berg 2008). The motion of a slender micro-swimmer through the fluid thus creates a hydrodynamic disturbance field similar to a force dipole on length scales much larger than the individual swimmer (Spagnolie & Lauga 2012). The cell structure, the number of flagella or cilia on the cell body and their arrangement collectively dictate the structure of the flow field around an individual swimmer.

The dipole hydrodynamic disturbance field is also observed in other bacteria such as Bacillus subtilis and in algae like Chlamydomonas reinhardtii. The far-field hydrodynamics of certain swimmers such as E. coli and B. subtilis involves drawing in fluid from the sides and expelling it along the longitudinal axis connecting the head and tail (figure 1a). They are hence classified as pusher-type swimmers due to the extensile nature of the resulting flow. In the case of C. reinhardtii, the synchronous beating of the twin flagella results in an inverse effect, wherein the fluid is drawn in parallel to the cell body and expelled along the sides of the swimmer (Guasto, Johnson & Gollub 2010; Garcia, Rafaï& Peyla 2013), as shown in figure 1(b). These swimmers are in turn termed as pullers, due to the contractile flow field that they generate through their motility. The difference in the sign of the dipole associated with the swimmer aids in classifying them as pusher- and puller-type swimmers. The magnitude of the dipole disturbance field created varies based on the cell structure and swimming speed of the organism.

Figure 1. A qualitative picture of the dipole hydrodynamic field of (a) pusher type swimmer such as E. coli or B. subtilis and (b) puller type swimmer such as the algae C. reinhardtii.

Swimmers tend to travel along fairly straight paths for long intervals of time, punctuated with short intervals involving a sudden change in their orientation. The sudden change is a consequence of one or more of the rotating flagella or cilia undergoing a switch in the direction of rotation. The resulting sharp change in orientation of the swimmer is due to the entire bundle of flagella or cilia coming undone and the swimmer itself coming to a momentary stop. This phenomenon is termed as tumbling (Berg 2008) and is a characteristic of some micro-swimmers which allows them to sample the surrounding fluid searching for sustenance. For long times this run–tumble behaviour of the bacterium is reminiscent of a random walk and gives rise to a diffusivity that is purely athermal in nature. Experiments by Rivero et al. (1989) and Chen, Ford & Cummings (2003) have revealed the preferential swimming of these organisms in a nutrient-rich medium, wherein the organism biases its random-walk by reducing its tumbling frequency when swimming along the attractant gradient. The work by Kalinin et al. (2009) has further revealed that the magnitude of the preferential swimming velocity is dependent on the gradient of the logarithmic attractant concentration.

The preferential swimming along the attractant gradient plays a crucial role in dictating the instability of an active suspension. For example, in the experiments by Dombrowski et al. (2004) and Tuval et al. (2005) involving bioconvection is a sessile drop, the aerotactic bacterium B. subtilis exhibits a preferential swimming towards the oxygen-rich interface. The observed large-scale coherence is attributed to the hydrodynamic interactions associated with the swimmer motility influenced by the interplay between aerotaxis and buoyancy. Tuval et al. (2005) have gone on to showcase the formation of a hydrodynamic vortex that aids in transporting oxygen from the interface into the drop through numerical simulations by modelling the interplay between aerotaxis and buoyancy effects using the Keller–Segel equations (Keller & Segel 1970). The simulations and experiments from Tuval et al. (2005) are of particular interest, as they yield valuable information pertaining to the orientation distribution of swimmers within the drop through the computed attractant field. Sokolov et al. (2009) performed experiments on thin liquid films containing suspensions of the aerotactic B. subtilis to determine a critical film thickness for the onset of collective motion. As a consequence of the preferential swimming, a clustering near the oxygen-rich interface resulting in the development of depletion layers is seen in both Tuval et al. (2005) and Sokolov et al. (2009).

Simha & Ramaswamy (2002) derived the hydrodynamic equations for a collection of self-propelled particles using general continuum principles and obtained an expression for the active stress – proportional to a dyadic product involving the orientational order parameter field. The continuum level modelling of an active suspension has been attempted by Saintillan & Shelley (2008) and Subramanian & Koch (2009) by using a Smoluchowski-type equation to describe the evolution of a probability distribution function $\varOmega (\boldsymbol {x},\boldsymbol {p},t)$, at a position $\boldsymbol {x}$ and orientation $\boldsymbol {p}$ of an active particle, associated with the instantaneous micro-structure of the suspension. The active stress in a suspension of swimmers was modelled by drawing motivation from the formulation given by Pedley & Kessler (1990), by building on the earlier studies on suspensions of passive anisotropic microstructures (Batchelor 1970; Hinch & Leal 1972). Subsequent works involving nonlinear simulations and stability calculations in the long-wave regime by Saintillan & Shelley (2013) have revealed the coupling between the active stresses within the suspension and the orientation field of the swimmer. Lushi, Goldstein & Shelley (2012) introduced an attractant within the system with its concentration field governed by an advection–diffusion equation. The linear stability of such a system was probed and the nonlinear evolution equations were solved with the inclusion of the effect of local gradients in attractant concentration on the orientation distribution of the swimmers. The more recent work by Lushi (2016) further accounts for the correlations associated with tumbling events of swimmers and probes the effects of anisotropy of the swimmer orientation field.

Motivated by the instability observed in the experiment performed by Sokolov et al. (2009) involving B. subtilis in a suspended thin film, Kasyap & Koch (2012) presented the formulation for an active stress-driven instability with an imposed attractant gradient. They showed through a long-wave stability calculation that the preferential swimming along the attractant gradient results in an anisotropy in the active stress that drives an instability in a suspension of pusher-type swimmers. The subsequent work by Kasyap & Koch (2014) further probes the stability characteristics in the finite wavenumber regime. We have recently extended the work by Kasyap & Koch (2014) by considering the effects of auto-chemotaxis Murugan & Roy (2022), wherein the evolution of the attractant field is dictated by effects of convection, diffusion and attractant secretion by the swimmers. We found that the self-signalling among swimmers through the secretion of an attractant is capable of generating an instability in a suspension of pullers, contrary to the calculation by Kasyap & Koch (2012, 2014) wherein it is shown that suspensions of pullers are unconditionally stable.

The ubiquitous nature of collective behaviour such as biofilm formation and swarming, along with their increased relevance near interfaces has drawn much attention to confined active matter suspensions. The work by Angelini et al. (2009) involving B. subtilis revealed the formation of surfactant waves which aid in the spreading of the bacterial film. Seminara et al. (2012) showed through experiments and theory that the presence of water-consuming extra-cellular matrix (ECM), typically seen in bio-films, creates an osmotic flux resulting in the inflow of water from the agar substrate which drives the film spreading. Both studies eliminate the possibility of the swimmer motility contributing to the film spreading. On the other hand, Sankararaman & Ramaswamy (2009) have formulated the flow dynamics for a thin film of polar self-propelled particles to showcase the instability arising as a consequence of the coupling between the active stress originating from the swimmer motility, the polarity and the free interface. Joanny & Ramaswamy (2012) described the spreading of an active droplet driven by stresses generated by the self-propelled particles contained within the drop. In both cases, the shape of the interface plays the role of setting up the inhomogeneities in the orientation field which results in the generation of an active stress that aids the droplet or film spreading. More recently, Alonso-Matilla & Saintillan (2019) showcased an instability associated with the interface in a polar active film containing a suspension of pushers.

The literature on swarming and biofilm modelling is quite extensive and various approaches have been used to model the physics of the system. The papers by Seminara et al. (2012) and Trinschek, John & Thiele (2016) depict the spreading of a film forced by osmotic gradients that develop due to the presence of the ECM seen in biofilms. Tam et al. (2019) proposed a model with contributions to the stress arising from the non-divergence free nature of the velocity field due to the effects of cell proliferation and death. Bees et al. (2000, 2002) modelled the spreading as a consequence of wetting, with the inclusion of van der Waals forces, which also mathematically serves to regularize the singularity that forms at a no-slip boundary. Fauvart et al. (2012) showed through both theory and experiments the swarming behaviour of Pseudomonas aeruginosa driven by surface forces. The recent work by Kotian et al. (2020) includes the effects of van der Waals forces, in addition to solving an evolution equation for the swimmer density field and accounting for its effects on the surfactant concentration field. Kotian et al. (2020) wrote down a free energy functional pertaining to the long-range van der Waals forces and the surface forces generated by the Marangoni stresses in the system. Subsequent minimization of the free energy functional, to generate the governing equations of the system, has also been done by Trinschek et al. (2017), with an added contribution to the free energy arising from the mixing of solute and solvent phases. Trinschek, John & Thiele (2018) has also explored the combined effects of both surface forces and wetting on the phenomenon of spreading. The recent study by Srinivasan, Kaplan & Mahadevan (2019) explores the physics of both swarming and biofilm propagation. They show that the swarming behaviour is prominently driven by capillarity and relies on wetting agents to generate an osmotic flux, while biofilm propagation relies on the osmotic flux generated by the ECM.

The current work aims to explore the coupled physics of the active stresses present in a suspension of swimmers and the effect of having a free interface through the use of a continuum formulation. By formulating the problem of a thin-film flow driven by the active stresses in the system, we are particularly interested in the evolution of the free interface in response to the active stress in the system. We predict that the anisotropy in the active stress, that results from the imposed attractant gradient, would create a jump in viscous stress at the interface that drives a localized shear flow. In addition, by drawing motivation from the rich literature on the long-wave nonlinear evolution of thin films (Oron & Rosenau 1994; Oron, Davis & Bankoff 1997; Craster & Matar 2009), using a long-wave lubrication theory we wish to write down a coupled set of nonlinear equations governing the film evolution driven by the active stresses in the system.

In the current study, in § 2 we start by formulating the thin-film problem with an imposed attractant gradient. The presence of the free interface results in symmetry breaking and leads to distinctly different stability characteristics for pushers and pullers. In § 3 we present a linear stability theory for the thin-film active suspensions, with the flow within the film being driven by a combination of the normal stress difference in the suspension, and the jump in normal and tangential stresses at the interface. We have added elaborate toy problems that discuss in detail the physics brought in by each of these stresses in Appendices B, C and D. In § 3.1 we characterize the instability of the coupled system in the long-wave regime analytically using a regular perturbation expansion. After validating these results numerically using a finite wavenumber numerical solver, in § 3.2 we probe the stability of the system in the finite wavenumber regime. Subsequently, in § 4 we formulate a long-wave theory for the film evolution. A coupled set of nonlinear equations governing the film height and the gap-integrated number density are solved numerically. Aside from validating the growth rates of the perturbations from the results obtained in the linear theory, we use the results from the numerical simulations to further explain the mechanism of the two modes of instability.

2. Problem formulation

The system under consideration involves a thin film of bacterial suspension subjected to an attractant gradient aligned along the $z$-direction, as shown in figure 2. At the continuum level, the film hydrodynamics is incompressible and governed by the Navier–Stokes equation

(2.1)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}$$

(2.2)$$\begin{gather}\rho \left( \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u} \right)={-} \boldsymbol{\nabla} p + \mu \nabla^{2} \boldsymbol{u} + \boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{T^{B}}.\end{gather}$$

Here, $\boldsymbol{u}$ is the velocity field, $\mu$ is the fluid viscosity and $\boldsymbol{T^{B}}$ represents the stress on the fluid arising as a consequence of the hydrodynamic field corresponding to the multitude of swimmers within the suspension. The active stress $\boldsymbol {T^B}$ can be written as an orientational average of the stress contribution from each individual swimmer (Subramanian & Koch 2009)

(2.3)\begin{equation} \boldsymbol{T^B}(\boldsymbol{x},t) ={-}C \mu \, U_s L^2 \int \varOmega(\boldsymbol{x},\boldsymbol{p},t) \left( \boldsymbol{p}\boldsymbol{p} - \frac{\boldsymbol{I}}{3} \right) {\rm d}\boldsymbol{p}, \end{equation}

where $C$ represents the non-dimensional dipole moment, $U_s$ is the swimming speed of the organism, $L$ is the length of the swimmer and $\varOmega (\boldsymbol {x},\boldsymbol {p},t)$ represents the probability distribution function associated with finding a swimmer at a location $\boldsymbol {x}$ having an orientation $\boldsymbol {p}$. The distribution $\varOmega (\boldsymbol {x},\boldsymbol {p},t)$ for perfectly isotropic tumbles is governed by a kinetic equation as given by Subramanian & Koch (2009)

(2.4)\begin{equation} \frac{\partial \varOmega}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} \left[ ( U_s \, \boldsymbol{p} + \boldsymbol{u} ) \right] + \boldsymbol{\nabla}\boldsymbol{\cdot} \left( \dot{\boldsymbol{p}} \, \varOmega \right) + \left( \frac{\varOmega}{\tau} - \frac{1}{4 {\rm \pi}} \int \frac{\varOmega}{\tau} \, {\rm d}\boldsymbol{p} \right) = 0,\end{equation}

where $\tau$ is the persistence time associated with the run and tumble motion of an individual swimmer. At length and time scales much larger than $U_s \tau$ and $\tau$, assuming uncorrelated tumbles, the kinetic equation reduces to a balance between swimmers tumbling out and into a given orientation

(2.5)\begin{equation} \left( \frac{\varOmega}{\tau} - \frac{1}{4 {\rm \pi}} \int \frac{\varOmega}{\tau} \,{\rm d}\boldsymbol{p} \right) = 0. \end{equation}

The above integro-differential equation has been solved by Subramanian, Koch & Fitzgibbon (2011) through the decomposition of the probability distribution function into a spatio-temporally varying number density field $n(\boldsymbol {x},t)$ and an orientation dependent $f(\boldsymbol {p})$, $\varOmega (\boldsymbol {x},\boldsymbol {p},t) = n(\boldsymbol {x},t) \hspace {0.05cm} f(\boldsymbol {p})$. The functional form of the biphasic tumbling frequency $(\tau ^{-1})$ as established by Chen et al. (2003) and Rivero et al. (1989) can be linearized in the limit of weak chemotaxis, characterized by $\xi = \chi \, U_s \, G \ll 1$. Here, $\chi$ is the sensitivity of the bacterium to the presence of an attractant gradient and $G$ indicates the strength of the attractant gradient. This linearized form of the biphasic tumbling frequency is used to arrive at the following expression for $f(\boldsymbol {p})$:

(2.6)\begin{equation} f(\boldsymbol{p}) = \begin{cases} \dfrac{1}{4 {\rm \pi}} \left[ 1 + \left( \boldsymbol{p} \boldsymbol{\cdot} \boldsymbol{g} - \dfrac{1}{4} \right) \xi \right] & \text{if} \ \boldsymbol{p} \boldsymbol{\cdot}\boldsymbol{g} > 0,\\[10pt] \dfrac{1}{4 {\rm \pi}} \left[ 1 - \dfrac{1}{4} \xi \right] & \text{if} \ \boldsymbol{p} \boldsymbol{\cdot} \boldsymbol{g} < 0. \end{cases} \end{equation}

Here, $\boldsymbol {g} = g \, \hat {z}$ and $g = \pm 1$ dictates the direction of the attractant gradient with respect to the interface and $\hat {z}$ is the unit vector along the vertical direction and points towards the interface. For films of sufficiently small thicknesses, the dynamics of the attractant transport is diffusion dominated and gives rise to a linear variation of the attractant across the film gap (as shown in Appendix A of Kasyap & Koch 2014). This gives rise to a constant value of the attractant gradient within the domain. The distribution function can be used to determine the active-stress tensor which is related to the single particle stresslet $(\boldsymbol {S})$ (Batchelor 1970)

(2.7)\begin{equation} \boldsymbol{T^{B}} = n \boldsymbol{S}, \quad \text{where}\ \boldsymbol{S} ={-}\frac{C}{16}\mu U_s L^{2} \xi \begin{bmatrix} - \dfrac{1}{3} & 0 & 0 \\ 0 & - \dfrac{1}{3} & 0 \\ 0 & 0 & \dfrac{2}{3} \end{bmatrix}. \end{equation}

Taking a zeroth moment of the kinetic equation (2.4), yields an advection–diffusion equation for the conservation of the density field $n(\boldsymbol {x},t)$

(2.8)\begin{equation} \frac{\partial n}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} \left[ ( \boldsymbol{g} \, U_{0} + \boldsymbol{u} ) n - D \, \boldsymbol{\nabla} n \right] = 0, \end{equation}

where $U_0 = \frac {1}{6} \, U_s \, \xi$ is the mean velocity of the swimmer along the attractant gradient and $D$ is the athermal diffusivity arising from the run and tumble motion of the swimmers. The velocity field pertaining to the fluid flow within the film is represented by $\boldsymbol {u} = (u,w)$ (see figure 2). Detailed derivations of the above continuum equations can be found in Kasyap & Koch (2014). The boundary conditions corresponding to the above system comprise of no-slip, no-penetration and no-flux conditions at the bottom wall ($z = 0$)

(2.9a)$$\begin{gather} u = 0, \end{gather}$$

(2.9b)$$\begin{gather}w = 0, \end{gather}$$

(2.9c)$$\begin{gather}D \frac{\partial n}{\partial z} - U_0 g n = 0. \end{gather}$$

The boundary conditions at $z = h$, involve the kinematic equation for the height evolution, the no-flux condition for the density field

(2.10a)$$\begin{gather} w = h_t + u h_x, \end{gather}$$

(2.10b)$$\begin{gather}D \left( \frac{\partial n}{\partial z} - h_x \frac{\partial n}{\partial x} \right) - U_0 g n = 0, \end{gather}$$

and the continuity of normal stress and tangential stress

(2.10c)\begin{align} &\left.\begin{gathered} \hat{\boldsymbol{n}} \boldsymbol{\cdot} \boldsymbol{T} \boldsymbol{\cdot} \hat{\boldsymbol{n}} ={-}\gamma \boldsymbol{\nabla}\boldsymbol{\cdot} \hat{\boldsymbol{n}}\\ \implies \, \left\{ p - p_{atm} \right\} -2 \mu \left\{ \frac{\partial w}{\partial z} \left( 1 - h_x^2 \right) - h_x \left( \frac{\partial u}{\partial z} + \frac{\partial w}{\partial x} \right) \right\} \frac{1}{1+h_x^2}\\ +\,\frac{\gamma h_{xx}}{(1 + h_x^2)^{3/2}} =\left\{ h_x^2 T_{xx}^B + T_{zz}^B \right\} \frac{1}{1+h_x^2}, \end{gathered}\right\} \end{align}

(2.10d)\begin{align} &\left.\begin{gathered} \hat{\boldsymbol{n}} \boldsymbol{\cdot} \boldsymbol{T} \boldsymbol{\cdot} \boldsymbol{\hat{t}} = 0\\ \implies \, \mu \left( \frac{\partial u}{\partial z} + \frac{\partial w}{\partial x} \right) \left( 1 - h_x^2 \right) + 2 \mu h_x \left( \frac{\partial w}{\partial z} - \frac{\partial u}{\partial x} \right) = \left( T_{xx}^B - T_{zz}^B \right) h_x, \end{gathered}\right\} \end{align}

where $\gamma$ is the surface tension at the film-air interface, $\hat{\boldsymbol{n}}$ and $\hat{\boldsymbol{t}}$ represent the unit normal and tangent vectors at the interface, and $\boldsymbol {T}$ represents the total stress in the suspension

(2.11a–c)\begin{align} \hat{\boldsymbol{n}} = \frac{1}{\sqrt{h_x^2+1}}\left[{-}h_x , 1 \right], \quad \boldsymbol{\hat{t}} = \frac{1}{\sqrt{h_x^2+1}}\left[ 1 , h_x \right], \quad \boldsymbol{T} ={-} \boldsymbol{\nabla} p + \mu \nabla^{2} \boldsymbol{u} + \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{T^{B}}. \end{align}

In the above equations, $h_t$ and $h_x$ depict the temporal and spatial derivative along the horizontal direction of the interface height. Equations (2.1), (2.2) and 2.8 are non-dimensionalized using the following scaling, with an asterisk * used to indicate the non-dimensional variables

(2.12a–f)\begin{equation} x^* = \frac{x}{H} \quad z^* = \frac{z}{H} \quad u^* = \frac{u}{U_0} \quad w^* = \frac{w}{U_0} \quad t^* = t \, \frac{D}{H^2} \quad p^* ={-} \frac{p}{\frac{3}{2} \, S_{zz} \langle n \rangle }, \end{equation}

and the non-dimensional numbers – Péclet ($Pe$), capillary ($Ca$), Schmidt ($Sc$) and Reynolds ($Re$) numbers

(2.13a–e)\begin{equation} Pe = \frac{U_0 H}{D} \quad \beta = \frac{3}{8}C \langle n \rangle L^2 H \quad Ca = \frac{\mu U_0}{\gamma} \quad Sc = \frac{\mu}{\rho D} \quad Re = \frac{\rho U_0 H }{\mu}, \end{equation}

with $\beta$ being an additional non-dimensional parameter that can be viewed as the ratio of bacterial stress to the viscous stress (Kasyap & Koch 2012). We would like to mention the presence of a typographical error in the definition of $\beta$ in Murugan & Roy (2022). Here, $\langle n \rangle$ is the average concentration of swimmers in the suspension. The equations can be written in the following non-dimensional form with the asterisk * being dropped for the sake of convenience

(2.14a)$$\begin{gather} \frac{\partial u}{\partial x} + \frac{\partial w}{\partial z} = 0, \end{gather}$$

(2.14b)$$\begin{gather}\frac{1}{Sc} \frac{\partial u}{\partial t} + Re \left( u \frac{\partial u}{\partial x} + w \frac{\partial u}{\partial z} \right)={-} \beta \frac{\partial p}{\partial x} + \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial z^2} + \frac{\beta}{3} \frac{\partial n}{\partial x} , \end{gather}$$

(2.14c)$$\begin{gather}\frac{1}{Sc} \frac{\partial w}{\partial t} + Re \left( u \frac{\partial w}{\partial x} + w \frac{\partial w}{\partial z} \right)={-} \beta \frac{\partial p}{\partial z} + \frac{\partial^2 w}{\partial x^2} + \frac{\partial^2 w}{\partial z^2} - \frac{2 \beta}{3} \frac{\partial n}{\partial z} , \end{gather}$$

(2.14d)$$\begin{gather}\frac{1}{Pe} \frac{\partial n}{\partial t} + \frac{\partial ( u n )}{\partial x} + \frac{\partial ( w n )}{\partial z} + g \frac{\partial n}{\partial z} - \frac{1}{Pe} \frac{\partial^2 n}{\partial x^2} - \frac{1}{Pe} \frac{\partial^2 n}{\partial z^2} = 0. \end{gather}$$

The value of surface tension at an air–water interface is $\gamma \approx 71.97 \ \textrm {mN} \ \textrm {m}^{-1}$. The mean velocity along the attractant-gradient is typically 10 % of the swimming speed, $U_0 \approx 3\ \mathrm {\mu } \textrm {m}\ \textrm {s}^{-1}$ for B. subtilis and $U_0 \approx {10}{-}{20} \ \mathrm {\mu } \textrm {m}\ \textrm {s}^{-1}$ for C. reinhardtii (Guasto et al. 2010). The motility of the micro-organisms in the suspension and the nature of the ambient medium where these organisms are present (such as the presence of biofilms) typically gives rise to values of viscosity in the range of $\mu \approx 10^{-2}\unicode{x2013}10^3 \ \textrm {Pa} \ \textrm {s}$ (Horvat et al. 2019). This yields a capillary number in range of $Ca \approx 10^{-7}\unicode{x2013}10^{-1}$. The typical diffusivity of a bacterial swimmer such as E. coli is $D\approx 10^{-10} \ \textrm {m}^2\ \textrm {s}^{-1}$ (Brenner, Levitov & Budrene 1998). Thus the Schmidt and the Reynolds numbers for the system turn out to be, $Sc= {O}(10^{4})$ and $Re= {O}(10^{-4})$, respectively, which allows us to neglect the inertial terms in the momentum equation. For a film thickness of $H \approx 100 \ \mathrm {\mu } \textrm {m}$ (as seen in Sokolov et al. 2009) and a mean swimming speed of $3 \ \mathrm {\mu } \textrm {m}\ \textrm {s}^{-1}$ (the mean swimming speed along the attractant is typically 10 % of the swimming speed), the Péclet number turns out to be $Pe \approx 3$. The puller-type swimmer C. reinhardtii has a swimming speed in the range $100\unicode{x2013}200 \ \mathrm {\mu } \textrm {m}\ \textrm {s}^{-1}$ (Guasto et al. 2010) and the corresponding value of $Pe$ for such a suspension would lie in the range $10\unicode{x2013}20$. The value of the dipole moment for an E. coli bacterium is $C\approx 0.57$, whereas $C\approx 4.11$ for C. reinhardtii (Subramanian & Koch 2009; Guasto et al. 2010). The length of the swimmer cells are $L\approx 12 \ \mathrm {\mu } \textrm {m}$ and $L\approx 19.5 \ \mathrm {\mu } \textrm {m}$. For an average concentration of swimmers, $\langle n \rangle \approx 2 \times 10^{10} \ \textrm {mm}^{-1}$ (Sokolov et al. 2009), the value of $\beta \approx 60$ for E. coli and $\beta \approx -1170$ for C. reinhardtii.

Figure 2. Schematic of the thin-film problem with an attractant gradient pointed towards the interface. The swimmers bias their random walk so as to align their mean swimming orientation along the gradient. The velocity field within the film is represented as $\boldsymbol {u} = (u,w)$. The normal and tangential unit vectors associated with the interface are represented by $\hat {\boldsymbol {n}}$ and $\hat {\boldsymbol {t}}$.

Experiments by Sokolov et al. (2009) involving an aerotactic suspension of bacteria in a quiescent liquid film revealed an instability, wherein they observed the emergence of large-scale convective motions beyond a critical film thickness. Kasyap & Koch (2012) modelled the observed hydrodynamic instability by considering an active stress driven flow within a channel, with the anisotropy in the stress field arising from the biased random walk of the individual swimmers. In the current work, we replace the top wall with a free interface. The presence of the free interface results in top–bottom symmetry breaking and leads to different stability characteristics with an inversion in the sign of the imposed attractant gradient. The non-dimensionalized boundary conditions at the bottom wall are

(2.15a)$$\begin{gather} u = 0, \end{gather}$$

(2.15b)$$\begin{gather}w = 0, \end{gather}$$

(2.15c)$$\begin{gather}\frac{\partial n}{\partial z} - Pe \, g n = 0. \end{gather}$$

Similarly, the non-dimensionalized boundary conditions at the interface ($z=h$) are

(2.16a)$$\begin{gather} w = \frac{1}{Pe} h_t + u h_x, \end{gather}$$

(2.16b)$$\begin{gather}\frac{\partial n}{\partial z} - Pe \, g n + h_x \frac{\partial n}{\partial x} = 0, \end{gather}$$

and the continuity of normal stress and tangential stress

(2.16c)\begin{align}& \left.\begin{gathered} \hat{\boldsymbol{n}} \boldsymbol{\cdot} \boldsymbol{T} \boldsymbol{\cdot} \hat{\boldsymbol{n}} ={-}\frac{\boldsymbol{\nabla} \boldsymbol{\cdot} \hat{\boldsymbol{n}}}{Ca}\\ \implies \, \beta \left\{ p - p_{atm} \right\} -2 \left\{ \frac{\partial w}{\partial z} \left( 1 - h_x^2 \right) - h_x \left( \frac{\partial u}{\partial z} + \frac{\partial w}{\partial x} \right) \right\} \frac{1}{1+h_x^2}\\ +\, \frac{h_{xx}}{Ca} \frac{1}{(1 + h_x^2)^{3/2}} = \left\{ h_x^2 T_{xx}^B + T_{zz}^B \right\} \frac{1}{1+h_x^2} = \frac{1}{3} \beta n (h_x^2 - 2)\frac{1}{1+h_x^2} , \end{gathered}\right\} \end{align}

(2.16d)\begin{align}& \left.\begin{gathered} \hat{\boldsymbol{n}} \boldsymbol{\cdot} \boldsymbol{T} \boldsymbol{\cdot} \boldsymbol{\hat{t}} = 0 \\ \implies \, \left( \frac{\partial u}{\partial z} + \frac{\partial w}{\partial x} \right) \left( 1 - h_x^2 \right) + 2 h_x \left( \frac{\partial w}{\partial z} - \frac{\partial u}{\partial x} \right) = \left( T_{xx}^B - T_{zz}^B \right) h_x = \beta n h_x. \end{gathered}\right\} \end{align}

The above stress boundary conditions highlight the role of active stresses at the interface. The viscous tangential stresses are not zero at the interface; they are proportional to the ‘active’ normal stress difference ($N_1^B=T_{xx}^B-T_{zz}^B$). This modification is reminiscent of the Marangoni stresses due to surface tension gradients (Leal 2007). One must note that suspensions of passive non-Brownian particles also exhibit normal stresses (Guazzelli & Pouliquen 2018). Contrary to the scenario presently under consideration, the normal stresses for the passive case would only exist in the presence of a background shear flow, depending on the local shear rate. Therefore for particle-laden free surface flows, normal stresses would vanish at the interface due to the stress-free boundary condition (Dhas & Roy 2022). Thus, microorganisms’ biased swimming due to an attractant leads to an anisotropic stress, consequently altering the tangential stress balance. We will show later, in Appendix C, that this modification to the tangential stress boundary condition is responsible for the new interfacial instability.

3. Linear stability analysis

To determine the stability characteristics of the system, we linearize the equations using a normal mode form

(3.1)\begin{equation} \left\{ u , w , p , n , h \right\} = \left\{ 0 , 0 , \bar{p}(z) , \bar{n}(z) , 1 \right\} + \left\{ \tilde{u}(z) , \tilde{w}(z) , \tilde{p}(z) , \tilde{n}(z) , \tilde{h} \right\} \exp({{\rm i}kx + \sigma t}). \end{equation}

As is typical in a linear theory, the perturbation amplitudes are assumed to be small compared with the base state variables of the system

(3.2a–e)\begin{equation} \tilde{n} \ll \bar{n} \quad \tilde{p} \ll \bar{p} \quad \tilde{u} \ll 1 \quad \tilde{w} \ll 1 \quad \tilde{h} \ll 1. \end{equation}

The base state is a quiescent suspension, with the number density distribution being governed by a balance between the flux associated with the preferential migration along the attractant gradient and the diffusive flux arising from unbiased random walk of the swimmers. The base state pressure field is thus a function of the local number density

(3.3a,b)\begin{equation} \bar{n} = \frac{g \, Pe \, e^{g \, Pe \, z}}{e^{g \, Pe } - 1}, \quad \bar{p} = p_{atm} - \frac{2}{3} \bar{n}. \end{equation}

The equations governing the perturbation are

(3.4)$$\begin{gather} {\rm i} k \tilde{u} + \mathcal{D} \,\tilde{w} = 0, \end{gather}$$

(3.5)$$\begin{gather}- \beta {\rm i} k \tilde{p} + \left\{ \mathcal{D}^2 - k^2 \right\} \tilde{u} + \frac{\beta}{3} {\rm i} k \tilde{n} = 0, \end{gather}$$

(3.6)$$\begin{gather}- \beta \mathcal{D} \, \tilde{p} + \left\{ \mathcal{D}^2 - k^2 \right\} \tilde{w} - \frac{2 \beta}{3} \mathcal{D} \, \tilde{n} = 0, \end{gather}$$

(3.7)$$\begin{gather}\frac{\sigma \tilde{n}}{Pe} + {\rm i} k \bar{n}\, \tilde{u} + \mathcal{D} ( \tilde{w} \bar{n} ) + g \mathcal{D} \, \tilde{n} - \frac{1}{Pe} \left\{ \mathcal{D}^2 - k^2 \right\} \tilde{n} = 0, \end{gather}$$

where $\mathcal {D} = {\textrm {d}}/{\textrm {d} z}$ represents the derivative along the vertical direction. Eliminating $\tilde {u}$ and $\tilde {p}$, the momentum equations reduce to a form that is reminiscent of the Orr–Sommerfeld equation, but devoid of inertia and being forced by an active stress arising from a gradient in the number density

(3.8a)$$\begin{gather} \left[ \mathcal{D}^4 - 2 k^2 \mathcal{D}^2 + k^4 \right] \tilde{w} + \beta k^2 \mathcal{D} \tilde{n} = 0, \end{gather}$$

(3.8b)$$\begin{gather}\left[ \mathcal{D}^2 - g \, Pe \, \mathcal{D} - \left( \sigma + k^2 \right) \right] \tilde{n} - g \, Pe^2 \bar{n} \tilde{w} = 0. \end{gather}$$

Similarly the linearized boundary conditions are

(3.9)\begin{equation} \left.\begin{aligned} \tilde{w} = 0 \\ \mathcal{D} \tilde{w} = 0 \\ \mathcal{D} \tilde{n} - g \, Pe \, \tilde{n} = 0 \end{aligned}\right\rbrace\text{ at}\ {z = 0} \qquad \left.\begin{aligned} \mathcal{D}^2 \,\tilde{w} + k^2 \tilde{w} = \beta k^2 \bar{n} \tilde{h} \\ \mathcal{D}^3 \,\tilde{w} - 3k^2 \mathcal{D} \,\tilde{w} + \beta k^2 \tilde{n} = \dfrac{k^4 \tilde{h}}{Ca} \\ \dfrac{\sigma \tilde{h}}{Pe} = \tilde{w} \\ \mathcal{D} \,\tilde{n} - g \, Pe \, \tilde{n} = 0 \end{aligned}\right\rbrace \text{ at}\ {z = 1.} \end{equation}

The perturbed fluid flow described by (3.8)–(3.9) is driven by a combination of the stresses generated by swimmer density and interface perturbations to the thin-film suspension. The swimmer density perturbations give rise to an anisotropic active stress by (3.8)–(3.9) throughout the entire suspension that results in the creation of a normal stress difference. The instabilities resulting from this mode can be studied by perturbing solely the swimmer density field as seen in Kasyap & Koch (2012). In Appendix B we detail the effects of a wall with slip on the stability calculation by Kasyap & Koch (2012), wherein we find that the hindrance to flow offered by a no-slip wall serves to dampen the instabilities when compared with a slippery wall.

The perturbations to the interface gives rise to a jump in the viscous stress across the interface resulting in a local shear flow. This mode of flow instability can be studied by perturbing solely the interface and not the swimmer density field, as shown in the toy problem in Appendix C. However, these two toy problems are only capable of destabilizing suspensions of pushers. The flow instability seen in suspensions of pullers (discussed in the following §§ 3.1 and 3.2) arises from the coupling between the swimmer density and interface perturbations, wherein the perturbations to the swimmer density field give rise to a jump in the normal stress across the interface (see the normal stress balance in (3.9)). This coupling is better analysed in detail through the toy problem formulated in Appendix D wherein we isolate the role of density perturbations in creating a jump in normal stress at the interface. In the following subsection we show that the signatures associated with these three toy problems that can be observed clearly through a long-wave stability calculation of the system.

3.1. Long-wave stability analysis of the system

In this section we will carry out a long-wave linear stability analysis of the system. By perturbing the system using long waves $(k \ll 1)$ it is possible to asymptotically solve the governing stability equations and extract some physical insights into the mechanism associated with the instabilities observed in the system. Thus, a perturbation expansion of the system variables is performed by treating $(k)$ as a small parameter

(3.10a)$$\begin{gather} \tilde{w} = \tilde{w}_0 + k^2 \tilde{w}_2 + k^4 \tilde{w}_4 + \dots, \end{gather}$$

(3.10b)$$\begin{gather}\tilde{n} = \tilde{n}_0 + k^2 \tilde{n}_2 + k^4 \tilde{n}_4 + \dots, \end{gather}$$

(3.10c)$$\begin{gather}\sigma = \sigma_0 + k^2 \sigma_2 + k^4 \sigma_4 + \dots . \end{gather}$$

The interface height $h$ has not been expanded as a series in $k^2$, due to the results of the stability calculation being invariant to such an expansion, as has been shown in Smith (1990). As (3.8) and (3.9) are comprised of only even powers of the wavenumber $k$, it is sufficient to consider only even powers in the perturbation expansion above. The equations at ${O}(1)$ are

(3.11a)$$\begin{gather} \mathcal{D}^4 \tilde{w}_0 = 0, \end{gather}$$

(3.11b)$$\begin{gather}\mathcal{D}^2 \tilde{n}_0 - g \, Pe \,\mathcal{D} \tilde{n}_0 - \sigma_0 \tilde{n}_0 - g \, Pe^2 \, \bar{n} w_0 = 0, \end{gather}$$

with boundary conditions

(3.11c)\begin{equation} \left.\begin{aligned} \tilde{w}_0= 0\\ \mathcal{D} \tilde{w}_0= 0\\ \mathcal{D} \tilde{n}_0 - g \, Pe \, \tilde{n}_0= 0 \end{aligned}\right\rbrace \text{ at}\ {z = 0} \qquad \left.\begin{aligned} \mathcal{D}^2 \tilde{w}_0&= 0\\ \mathcal{D}^3 \tilde{w}_0&= 0\\ \frac{\sigma_0 \tilde{h}}{Pe}&= \tilde{w}_0\\ \mathcal{D} \tilde{n}_0 - g \, Pe \, \tilde{n}_0&= 0 \end{aligned}\right\rbrace \text{ at}\ {z = 1}. \end{equation}

The solution reveals the absence of flow at ${O}(1)$ of the perturbation expansion with an exponentially varying density profile

(3.12a)\begin{equation} \tilde{w}_0 = 0, \quad \tilde{n}_0 = \tilde{m}_0 \frac{g \, Pe \, e^{g \, Pe \, z}}{e^{g \, Pe} - 1 }, \quad \sigma_0 = 0, \end{equation}

where $\tilde {m}_0$ is the gap-integrated number density $n_0$. The equations at ${O}(k^2)$ are

(3.13a)$$\begin{gather} \mathcal{D}^4 \tilde{w}_2 + \beta \mathcal{D} \tilde{w}_0 = 0, \end{gather}$$

(3.13b)$$\begin{gather}\mathcal{D}^2 \tilde{w}_2 - g \, Pe \, \mathcal{D} \tilde{w}_2 - \sigma_2 \tilde{w}_0 - \tilde{w}_0 - g \, Pe^2 \, \bar{n} \tilde{w}_2 = 0, \end{gather}$$

subject to the boundary conditions

(3.14)\begin{equation} \left.\begin{aligned} \tilde{w}_2&= 0\\ \mathcal{D} \tilde{w}_2&= 0\\ \mathcal{D} \tilde{n}_2 - g \, Pe \, \tilde{n}_2 &= 0 \end{aligned}\right\rbrace \text{ at}\ {z = 0}\qquad \left.\begin{aligned} \mathcal{D}^2 \tilde{w}_2 &= \beta \bar{n} \tilde{h}\\ \mathcal{D}^3 \tilde{w}_2 + \beta \tilde{n}_0 &= \tilde{h} / \widehat{Ca}\\ \frac{\sigma_2 \tilde{h}}{Pe} &= \tilde{w}_2 \\ \mathcal{D} \tilde{n}_2 - g \, Pe \, \tilde{n}_2 &= 0 \end{aligned}\right\rbrace \text{ at}\ {z = 1}. \end{equation}

A modified capillary number, $\widehat {Ca} = Ca / k^2$, is defined to be an ${O}(1)$ quantity. The analytical solution of the above system yields a quadratic equation governing the growth rate $\sigma$. We relegate the detailed expressions of velocity, number density and growth rate to Appendix A. The stability characteristics reveal a new mode of instability associated with the interface and a modification to the previously characterized density mode by Kasyap & Koch (2012). The positive and negative values of $\beta$ correspond to the pusher- and puller-type swimmers that have been approximated as extensile and contractile force dipoles respectively. Pullers do not trigger an instability through the density mode detailed by Kasyap & Koch (2012) or by the decoupled interface mode (no number density perturbations) discussed in Appendix C. For the coupled system, we find a new result, a suspension of pullers is unstable to long-wave perturbations at finite and large $Pe$ for $g=+1$.

The neutral curves in the second quadrant of figure 3 pertaining to a suspension of pushers $(\beta > 0)$ with $g=-1$, displays a switch in the mechanism of instability from an interface mode to a density mode with increasing $Pe$. The system is strongly coupled and ascertaining the mechanism of the instability of an individual mode is done by drawing a comparison with the velocity profiles obtained in the toy problems discussed in Appendixes B and C. The velocity field of a density mode displays a flip in the sign of the horizontal velocity $u$ across the film height, whereas the interface mode velocity profile is reminiscent of a flow generated by applying a shear stress at the free interface. The transition from the interface mode to the density mode is thus observed by tracking the velocity field associated with the unstable mode as shown in figure 4. At low values of $Pe$, the uniform distribution of the density field creates a jump in viscous stress that drives the interface mode of instability. The solution of the quadratic equation governing the growth rate, for the asymptote of $Pe \ll 1$ yields

(3.15)$$\begin{gather} \sigma^{\textit{Interface Mode}} =\left\{ \frac{1}{2} \beta - \frac{1}{3 \, \widehat{Ca}} \right\} k^2 Pe + \left\{ \frac{1}{4} g \beta \right\} k^2 Pe^2 + {O}(Pe^3), \end{gather}$$

(3.16)$$\begin{gather}\sigma^{\textit{Density Mode}} ={-}k^2 - \left\{ \frac{1}{8} g \beta \right\} k^2 Pe^2 + {O}(Pe^3). \end{gather}$$

At leading order, this yields a critical activity of $\beta = 2 / ( 3 \, \widehat {Ca} )$ for the onset of an instability. At large $Pe$, for $g=-1$, a large number of bacteria succeed in accumulating near the bottom wall due to the weak diffusive flux in the system. This greatly diminishes the strength of the interfacial stress which drives the interface mode. Thus, the instability of the system is plausible only through a density mode, with the active stress generated by the pusher-type swimmers reinforcing the density perturbation.

Figure 3. Regions of stability from a long-wave theory for a suspension of pushers ($\beta >0$) and pullers (${\beta <0}$). The first and fourth quadrants pertain to a system with an attractant gradient pointed towards the interface $(g=+1)$, while the second and third quadrant correspond to the system with the gradient pointed towards the bottom wall $(g=-1)$ for (a) $\widehat {Ca}=0.1$, (b) $\widehat {Ca}=1$ and (c) $\widehat {Ca}=10$.

Figure 4. Velocity fields for $\widehat {Ca} = 1$. Points A, B, C and D correspond to pusher-type swimmers with $\beta = 10$ and $Pe = 10$, 3.5, 1.0 and $5$ respectively. Point E corresponds to Puller-type swimmers with $\beta = -10$ and $Pe = 5$. In the plots of the velocity fields pertaining to the two modes in the problem, the red and green curves depict the velocity field associated with unstable and stable modes, respectively.

The first quadrant in figure 3 corresponds to a suspension of pushers that exhibit preferential swimming towards the interface $(g=+1)$. This leads to the presence of a finite number of bacteria near the interface regardless of the strength of diffusion flux in the system. At small $Pe$ the interfacial stress is strong enough to drive an instability. The critical activity for the instability is independent of the direction of the attractant gradient ($g$) for $Pe \ll 1$, as is seen from (3.15). With increasing $Pe$, the strength of the interfacial stress is amplified due to the increasing number of swimmers near the interface. In addition, the active stress associated with the pusher field becomes significant with increasing $Pe$, as the weak diffusive flux fails to dampen the density field perturbation.

The fourth quadrant in figure 3 corresponds to a suspension of pullers with $g=+1$. The contractile nature of the puller velocity field makes it stable to perturbations associated with the density field as seen from Kasyap & Koch (2012, 2014). It can also be seen from Appendix C (C3) that a suspension of pullers $(\beta <0)$ is stable to perturbations associated with the interface. In the coupled system, the hydrodynamic field arising from the active stress strives to reinforce the interface perturbations, whereas the interfacial stress in the case of pullers drives a flow that strengthens the density perturbations, rendering the coupled system unstable. This mechanism associated with the instability is further verified when looking at the nonlinear long-wave simulations for the film evolution in § 4.2.

The system consisting of a suspension of pullers with $g=-1$ is unconditionally stable as seen from the third quadrant of figure 3. The coupled dynamics of both the interfacial and active stress is required for the instability of a suspension of pullers. In this system configuration, the active stress is severely weakened in the low $Pe$ regime of strong diffusion, wherein gradients in the density field are strongly suppressed. The finite and large $Pe$ regime involves the migration of swimmers towards the attractant-rich bottom wall, thus weakening the interfacial stress which is dependent on the density field of swimmers at the interface.

3.2. Finite wavenumber stability characteristics

In the previous sections, we discussed the instability mechanisms of the interface and density modes using a long-wave theory. We will now probe the stability of a chemotactic thin film to disturbances of arbitrary wavelengths, analysing the modification of both the instabilities. Kasyap & Koch (2014) used a spectral method to numerically probe the stability characteristics at finite wavenumber ($k$), for the active suspension system in a channel. They employed basis functions originating from the homogeneous solution of the equations in (3.8) subject to boundary conditions associated with the channel. In the current problem, the free interface makes the calculation of a similar set of basis functions tedious. We employ Lagrange polynomials with the domain being discretized using the Chebyshev grid points in the same manner as detailed by Trefethen (2000) to probe the finite $k$ stability characteristics

(3.17a)\begin{equation} \tilde{w} = \sum_{j=0}^{N-1} L_{ij} \tilde{w}_j, \quad \tilde{n} = \sum_{j=0}^{N-1} L_{ij} \tilde{n}_j, \end{equation}

where $w_j$ and $n_j$ are the values of $w$ and $n$ evaluated at the Chebyshev grid points, $z_j = \cos (j {\rm \pi}/ N)$. The spectral solver is validated with the long-wave stability calculation as shown in figure 5. The solver is then used to probe the stability characteristics at finite wavenumber.

Figure 5. Validation of long-wave theory with the numerical calculation for $Pe=5$, $\beta =100$, $Ca=10^{-4}$, $k=0.01$.

3.2.1. Pushers ($\beta >0$) with $g=-1$

The long-wave stability calculation for pushers with $g=-1$ (figure 4) revealed the presence of interface and density modes of instability. The interface mode of instability is stabilized at high values of $k$ as stronger gradients in the interface height end up amplifying the effects of surface tension. This can be seen from figure 6 which shows the switch in mechanism associated with the most unstable mode that occurs with increasing $k$. Figure 7 depicts the change in regions of stability in the $k\unicode{x2013}\beta$ plane for various values of $Pe$.

Figure 6. For a suspension of pushers with $g = -1$, $Pe = 4$, $\beta = 100$ and $Ca = 10^{-3}$, figure (a) depicts the variation of the growth rate ($\sigma$) with the wavenumber ($k$), for the two unstable modes in the problem (red and blue lines) with an asymptote (dashed black line) providing a comparison with the long-wave theory discussed in § 3.1. Panels (b,c) depict the velocity field for different values of the wavenumber ($k$) pertaining to the growth rates of the two modes displayed in (a), with the in-line numbers depicting the appropriate wavenumber $(k)$ at which the eigenfunction is generated. In (b,c) the dashed curves indicate that the mode has stabilized at the corresponding inline wavenumber.

Figure 7. Neutral curves in the $k\unicode{x2013}\beta$ plane for pushers with $g=-1$ and $Ca=10^{-3}$ for (a) $Pe=0.1$, (b) $Pe=1.0$, (c) $Pe=4.0$ and (d) $Pe=10.0$.

At low values of $Pe$ a high value of activity ($\beta$) is required to trigger an instability through the density mode (already established in Kasyap & Koch (2012) and Appendix B). This is due to the strong diffusive flux in the system, which disallows gradients in the density field as seen in the long-wave calculation. This statement can be inferred by comparing figures 7(a) and 7(d) which correspond to $Pe=0.1$ and $Pe=10.0$, and showcases the absence and presence of finite wavenumber instabilities, respectively. It is more readily apparent when looking at the figure 8, wherein the neutral regions are plotted in the $Pe\unicode{x2013}\beta$ plane.

Figure 8. Neutral curves in the $Pe\unicode{x2013}\beta$ plane for pushers with $g=-1$ and $Ca=10^{-3}$ for (a) $k=0.04$, (b) $k=0.4$, (c) $k=1.0$ and (d) $k=8.0$.

The instability of the system in the small $Pe$ regime is thus primarily dominated by the interface mode which is particularly strong at low wavenumbers wherein the effects of surface tension are diminished. The subsequent transition wherein instabilities appear in the finite wavenumber regime as seen in figures 7(b)–7(c) is due to the increased prominence of the density mode, especially at high values of $k$ wherein the stronger gradients in the density field amplify the pusher velocity field and enhance the density perturbations. Figure 8 displays the regions of stability in the $Pe\unicode{x2013}\beta$ plane. In particular, the lack of instabilities in the large $k$ regime, for small and finite $Pe$ numbers be observed from figure 8(d).

3.2.2. Pushers ($\beta >0$) with $g=+1$

The attractant gradient pointing towards the interface causes a finite number of bacteria to be present at the interface regardless of the strength of diffusion in the system. The interface mode is driven by the jump in viscous stress and is proportional to the density of swimmers at the interface. At small $Pe$ the regime is diffusion dominated and the effect of preferential swimming arising from the direction of the attractant gradient is negligible. The interfacial stress varies as ${\sim}O(1)$. The neutral curves corresponding to low $Pe$ in figure 9 show that the critical activity $\beta$ required for the onset of instability is lower in the long-wave regime. This is consistent with the characteristics of the interface mode which stabilizes at finite wavenumber due to surface tension effects. At moderate and large $Pe$, the diminished effect of diffusion allows for a large number of swimmers to accumulate near the interface. The interfacial stress varies as ${\sim}O(Pe)$ for large $Pe$ which explains the reduction in the critical value of $\beta$ for the onset of instability as seen in 9. Furthermore, gradients in the density field become prominent owing to weak diffusion thus allowing the active stress to further attenuate the interface perturbations. Figure 10 reinforces this result as it is clearly seen from the $Pe\unicode{x2013}\beta$ plane that the large $Pe$ limit becomes unstable at much lower values of activity regardless of the value of $k$. The upward shift of the neutral curves with increasing $k$ is indicative of the increasing strength of surface tension at finite wavenumber.

Figure 9. Neutral curves in the ${k}{-}{\beta}$ plane for pushers with $(g=+1)$ and $Ca=10^{-3}$ for (a) $Pe=0.1$, (b) $Pe=2.0$ and (c) $Pe=6.0$.

Figure 10. Neutral curves in the $Pe\unicode{x2013}\beta$ plane for pushers with $(g=+1)$ and $Ca=10^{-3}$ for (a) $k=0.1$, (b) $k=3.5$ and (c) $k=8.0$.

3.2.3. Pullers $(\beta <0)$ with $g=+1$

Through the long-wave stability calculation in § 3.1 and the toy problem formulated in Appendix D, it was shown that the coupling between the perturbations applied to the swimmer density field and the interface height is crucial in destabilizing suspensions of pullers. The stability calculation for suspensions of pullers involving finite wavenumber perturbations reveal a similar trend, wherein the dampening of perturbations associated with either the swimmer density field or the interface height is sufficient to stabilize the system. This result can be seen observed from figure 11, wherein the regions of stability in the $Pe\unicode{x2013}\beta$ across different wavenumbers ($k$) display the universal trend of stabilizing with decreasing $Pe$. This is contrary to the behaviour of suspensions of pushers, wherein the lowering of the strength of diffusive flux in the system, merely serves to at best switch the mode of instability.

Figure 11. Neutral curves in the ${Pe}{-}{\beta}$ plane for pullers with $g=+1$ and $Ca = 10^{-3}$ for (a) $k=0.04$, (b) $k=0.4$ and (c) $k=2.0$.

A similar phenomenon occurs at large wavenumbers wherein the effects of surface tension are strong enough to dampen the perturbations to the interface, consequently stabilizing the suspension as a whole, as shown in figure 12. Figures 11 and 12 together showcase that the coupling between density and interface perturbations is the driving mechanism for the instability even at finite wavenumbers.

Figure 12. Neutral curves in the ${k}{-}{\beta}$ plane for pullers with $g=+1$ and $Ca = 10^{-3}$ for (a) $Pe=0.5$, (b) $Pe=2.0$ and (c) $Pe=6.0$.

To further understand the mechanism associated with the instability, it is necessary to look at the evolution of the system beyond the purview of a linear theory. In the next section, a set of nonlinear evolution equations for the system in the long-wave limit are formulated and solved.

4. Nonlinear long-wave theory of the film

The problem formulation in §§ 2 and 3 made use of a single length scale, the film height $H$, to scale lengths along the $x$ and $z$ directions. The essence of developing a long-wave theory for the film evolution relies on the disparity in length scales in the two directions – film height $H$ in the $z$-direction and wavelength of disturbances $\lambda$ in the $x$-direction. This is done through the definition of a film aspect ratio that is assumed asymptotically small $(\epsilon = H / L \ll 1)$. Thus, all lengths along the $z$ and $x$ direction are scaled with the film height $H$ and $H / \epsilon$, respectively. This approximation makes the velocities along the $x$ and $z$ directions differ by ${O}(\epsilon )$.

The $z$-momentum equation involves a balance between the hydrostatic pressure and the active stress, through which we obtain a scaling for the pressure field. The $x$-momentum equation involves a balance between the hydrostatic pressure, the active stress, and the viscous stress. The time scaling is chosen such that the unsteady term from the number density conservation equation balances the difference between the diffusive and advective flux of the density field

(4.1a–e)\begin{equation} x^* = x H \quad z^* = z \epsilon H \quad u^* = \epsilon U_0 u \quad w^* = \epsilon^2 U_0 w \quad t^* = t \epsilon^{{-}2} \frac{H^2}{D}. \end{equation}

With the long-wave scaling arguments presented above the number density and momentum conservation equations transform as

(4.2a)$$\begin{gather} \left[\frac{1}{Pe} \frac{\partial^2 n}{\partial {z}^{2}} - g \frac{\partial n}{\partial z} \right] + \epsilon^2 \left[ \frac{1}{Pe} \frac{\partial^2 n}{\partial {x}^2} - \frac{1}{Pe}\frac{\partial n}{\partial t} - \frac{\partial \left( u n \right) }{\partial x} - \frac{\partial \left( w n \right) }{\partial z} \right] = 0, \end{gather}$$

(4.2b)$$\begin{gather}\epsilon \left[ - \beta \frac{\partial p}{\partial x} + \frac{\partial^2 u }{\partial {z}^{2}} + \frac{\beta}{3} \frac{\partial n}{\partial {x}} \right] + \epsilon^3 \left[ \frac{\partial^2 u}{\partial {x}^{2}} \right] = 0, \end{gather}$$

(4.2c)$$\begin{gather}\left[ \frac{\partial p}{\partial z} + \frac{2}{3} \frac{\partial n }{\partial z} \right] + \epsilon^2 \left[ \frac{1}{\beta} \frac{\partial^2 w}{\partial {z}^{2}} \right] + \epsilon^4 \left[ \frac{1}{\beta} \frac{\partial^2 w}{\partial {x}^{2}} \right] = 0. \end{gather}$$

The boundary conditions at the bottom wall remain the same as in (2.15). At the interface, $z = h$, applying the scaling proposed by the long-wave theory in (4.1a–e), the boundary conditions in (2.16) transform as

(4.3a)$$\begin{gather} \left[\beta p + \frac{2}{3} \beta n + \frac{h_{xx}}{\widehat{Ca}} \right] + \epsilon^2 \left[ \beta p h_x^2 - \frac{1}{3} \beta n h_x^2 + 2 h_x \frac{\partial u}{\partial z} + \frac{h_{xx}}{\widetilde{Ca}} h_x^2 \right] = 0, \end{gather}$$

(4.3b)$$\begin{gather}\epsilon \left[ \frac{\partial u}{\partial z} - \beta n h_x \right] - \epsilon^2 \left[ h_x \frac{\partial u}{\partial z} \right] + \epsilon^3 \left[ \frac{\partial w}{\partial x} - h_x^2 \frac{\partial u}{\partial z} + h_x \frac{\partial w}{\partial z} \right] = 0, \end{gather}$$

(4.3c)$$\begin{gather}w = \frac{1}{Pe} h_t + u h_x, \end{gather}$$

(4.3d)$$\begin{gather}\left[\frac{\partial n}{\partial z} - Pe \, g n \right] + \epsilon^2 \left[ h_x \frac{\partial n}{\partial x} \right] = 0. \end{gather}$$

In the above set of equations, the subscripts $t$ and $x$ indicate the derivatives of the variable with respect to time ($t$) and the horizontal coordinate ($x$), respectively, with $\widetilde {Ca} = Ca / \epsilon ^2 = \mu U_0 / \epsilon ^2 \gamma$ denoting a modified capillary number. For a film aspect ratio $\epsilon = H / L \leq 0.1$, the maximum value of the modified capillary number $\widetilde {Ca} \approx 1$, and is assumed to be an ${O}(1)$ quantity.

4.1. Nonlinear evolution equations

Equation (4.2a) at ${O}(1)$ subject to the above no-flux conditions gives the form of the number density

(4.4)\begin{equation} n(x,z,t) = m(x,t) \frac{g \, Pe \, e^{g \, Pe \, z}}{ e^{g \, Pe \, h} - 1} = C(x,t) e^{g \, Pe \, z}, \end{equation}

where $m = \langle n(x,t) \rangle$ is the gap-integrated number density and the variable $C(x,t) = m ( x , t ) g \, Pe / ( e^{g \, Pe \, h} - 1 )$ is defined for the sake of algebraic simplicity. Here, $\langle \, \rangle = \int _{0}^{h} \,\textrm {d}z$, denotes a gap integration. Integrating the $z$-momentum equation (4.2b) at ${O}(1)$ yields an expression for the pressure, with the constant of integration being evaluated from the normal stress balance at the interface (4.3a)

(4.5)\begin{equation} p(z,x,t) ={-} \frac{h_{xx}}{\beta \widetilde{Ca}} - \frac{2 }{3}n(z,x,t) + p_{atm}. \end{equation}

Using the pressure field calculated above, the $x$-momentum equation (4.2b) at ${O}(\epsilon )$ yields a second-order differential equation governing $u$

(4.6)\begin{equation} \frac{\partial^2 u }{\partial {z}^{2}} = \left[ - \frac{h_{xxx}}{\widetilde{Ca}} - \beta \frac{\partial n}{\partial {x}} \right]. \end{equation}

The above equation can be integrated twice with the constants of integration being determined from the no-slip boundary condition at the bottom wall ($z=0$) and (4.3b), to obtain an expression for the horizontal velocity field $u$

(4.7)\begin{equation} u ={-} \frac{h_{xxx}}{\widetilde{Ca}} \left[ \frac{z^2}{2} - h z \right] - \frac{\beta C_x}{g \, Pe} \left[ \frac{e^{g \, Pe \, z} - 1 }{g \, Pe} - z e^{g \, Pe \, h} \right] + \beta h_x z C e^{g \, Pe \, h}. \end{equation}

The first term in the velocity field captures the effect of surface tension, which tries to suppress the perturbations to the interface. The second term arises from the pusher or puller stress on the fluid and is responsible for the density mode observed in the long-wave stability calculation. The third term is a consequence of the jump in viscous stress at the interface and is reminiscent of the Marangoni stress seen in Oron & Rosenau (1994). Having arrived at an expression for the horizontal velocity field, a gap integration performed on the number density equation (4.2a) at ${O}(\epsilon ^2)$ yields

(4.8)$$\begin{gather} m_t + \frac{\partial }{\partial x} \left\{ \int_0^h \left( Pe \, u n - \frac{\partial n}{\partial {x}} \right) {\rm d}z \right\} = 0, \end{gather}$$

(4.9)$$\begin{gather}h_t + \frac{\partial}{\partial x} \left\{\int_0^{h} Pe \, u \, {\rm d}z \right\} = 0. \end{gather}$$

The coupled set of equations are expanded and written using the expressions for $n$ and $u$ obtained in (4.4) and (4.7)

(4.10)\begin{align} &\frac{1}{Pe} h_t + \frac{\partial }{\partial x} \left[ \frac{h_{xxx} h^3}{3 \widetilde{Ca}} + \frac{1}{2} \beta h_x h^2 C e^{g \, Pe \, h} \right.\nonumber\\ &\quad \left.+ \,\frac{ \beta C_x}{2 g^3 \, Pe^3} \left( 2 + 2 g \, Pe \, h + e^{g \, Pe \, h} \left\{ g^2 \, Pe^2 \, h^2 - 2 \right\} \right) \right] = 0, \end{align}

(4.11)\begin{align} &\frac{1}{Pe} m_t + \frac{\partial }{\partial x} \left[ \frac{C h_{xxx} }{2 \, \widetilde{Ca} \, g^3 \, Pe^3} \left\{ 2 + 2 g \, Pe \, h + e^{g \, Pe \, h} \left( g^2 \, Pe^2 \, h^2 - 2 \right) \right\}\right.\nonumber\\ &\quad +\, \frac{\beta C C_x }{2 g^3 \, Pe^3} \left\{ 4 e^{g \, Pe \, h} - 1 + e^{2 g \, Pe \, h} \left( 2 g \, Pe \, h - 3 \right) \right\}\nonumber\\ &\quad \left. + \frac{ \beta h_x C^2 e^{g \, Pe \, h}}{g^2 \, Pe^2} \left\{ 1 + e^{g \, Pe \, h} \left( g \, Pe \, h - 1 \right) \right\} - C_x \frac{e^{g \, Pe \, h} - 1}{g \, Pe^2 } \right] = 0. \end{align}

4.2. Numerical solution of nonlinear long-wave equations

The coupled long-wave nonlinear system of equations (4.10) and (4.11) are solved numerically using a fourth-order Runge–Kutta time integrator with the spatial derivatives being computed using a second-order finite-difference scheme. Of the four possible combinations involving suspensions of pushers and pullers subjected to an attractant gradient along the negative or positive vertical axis, there are three configurations that exhibit long-wave hydrodynamic instabilities according to the linear theory derived in the previous section. We investigate the evolution of the unstable system through both mechanisms of instability unearthed in the long-wave linear stability calculation. The system of pushers subjected to a negative attractant gradient is of particular interest, since there is a switch in the mechanism of instability of the most unstable mode in the system with increasing $Pe$, as shown in figure 13.

Figure 13. (a) The regions of stability and (b) the variation of the growth rate with $Pe$ for $\beta = 5$ for a suspension of pushers with a negative attractant gradient at $\widehat {Ca} = 0.1$ and $k = 0.1$.

4.2.1. Pushers with a negative gradient

The long-wave calculation for a suspension of pushers subjected to an attractant gradient pointed towards the bottom wall reveals three distinct regions of instability in the $Pe\unicode{x2013}\beta$ plane as seen in figure 13. The intermediate region coloured in orange indicates the parameter space where both the modes are unstable. The nonlinear long-wave equations are solved numerically at the three points A, B and C depicted in figure 13 and the resulting time evolution of the interface height and gap-integrated density field is plotted in figure 14.

Figure 14. Nonlinear evolution of interface and pusher density field for $g=-1$, $\beta = 5$, $\widetilde {Ca} = 0.1$ and $k = 0.1$. Solid lines indicate the interface height $h(x,t)$ and dashed lines the gap-integrated density field $m(x,t)$ for (a) $Pe = 1.0$ (interface mode), (b) $Pe = 3.0$ (interface mode), (c) $Pe = 5.0$ (density mode), (d) $Pe = 5.0$ (interface mode) at different times as indicated in the curve inset.

The instability of the interface mode is driven by the jump in viscous stress which has been shown to be proportional to the swimmer concentration at the interface. At low values of $Pe$, the dominant diffusive flux in the number density conservation works to set up a uniform distribution of swimmers along the film height and strives to negate the biased migration of swimmers to accumulate near the bottom wall. This ensures the presence of a finite number of swimmers at the interface resulting in the generation of the required interfacial stress which favours the growth of the perturbation as seen in figures 14(a) and 14(b).

In the limit of large $Pe$, the negligible diffusive flux of the swimmer density field leads to a large concentration of swimmers accumulating near the bottom wall. This diminishes the jump in viscous shear stress at the interface and greatly reduces the growth rate associated with the interface mode of instability. On the other hand, the lack of diffusive flux fails to disperse the growing concentration of swimmers driven by the hydrodynamic velocity field of a pusher. Thus, we observe an instability arising from a mechanism similar to that detailed by Kasyap & Koch (2012), which we have termed as the density mode (figure 14c).

The streamlines in figure 15 depict the flow draining out from under the trough, for the case of the density mode (figure 14c). Furthermore, it is observed that the effect of surface tension works to flatten the perturbation to the interface but the more significant contributions arising from the active and interfacial stress succeed in deepening the trough in a manner that is reminiscent of a film rupture. The numerical simulations are restricted by the accuracy of the solver in resolving the field variables and their higher-order derivatives at the deepening trough of the interface.

Figure 15. Streamlines for the pusher density mode for an imposed attractant gradient along the negative $z$-direction $(g=-1)$ for $\beta = 5$, $\widetilde {Ca} = 0.1$ and $k = 0.1$ at time $t=140.97$.

4.2.2. Pushers with a positive gradient

By flipping the direction of the attractant gradient for a suspension of pushers, we expect a similar instability at low $Pe$ since the density profile tends towards being distributed uniformly across the film depth. Figures 16(a) and 16(b) display the effects of surface tension in suppressing the waves that form on the interface. Figure 16(c) depicts the saturation of the growth rate associated with the perturbation. It is interesting to note that the density field closely shadows the interface evolution. This further supports the argument that the strong effect of diffusion reduces the effect of clustering of the swimmers near an attractant-rich boundary and that the variations in the gap-integrated density field are solely caused by the deformation of the interface.

Figure 16. Nonlinear evolution of interface and pusher density field for an imposed attractant gradient along the positive $z$-direction. Solid lines indicate the interface height $h(x,t)$ and dashed lines the gap-integrated density field $m(x,t)$ for $Pe = 0.1$, $\beta = 10$, $k = 0.1$ and (a) $\widetilde {Ca} = 0.01$, (b) $\widetilde {Ca} = 0.001$ at different times as indicated in the figure inset. Panel (c) depicts the growth rate of the perturbations with time.

At finite values of $Pe$, we have an increasing concentration of swimmers accumulating near the attractant-rich interface, which serves to strengthen the interfacial stress. The evolution of the system for finite $Pe$ is depicted in figures 17(a) and 17(b). The perturbation grows unboundedly (figure 17c) and leads to the divergence of the numerical solution. The weak diffusive flux in the system fails to suppress the growth of the density perturbation. Furthermore, the clustering of swimmers near the interface further amplifies the interface perturbation. These effects can be seen by observing the flow streamlines and velocity contours displayed in figure 18 which shows that the active-stress-driven flow and the interfacial stress serve to amplify the interface perturbation. It is for this reason that we do not observe a saturation of the perturbation growth as seen in the low $Pe$ regime. The effect of surface tension merely dampens the growth rate but ultimately fails to saturate the growth of the perturbation.

Figure 17. Nonlinear evolution of interface and pusher density field for an imposed attractant gradient along the positive $z$-direction. Solid lines indicate the interface height $h(x,t)$ and dashed lines the gap-integrated density field $m(x,t)$ for $Pe = 5$, $\beta = 2$, $k = 0.1$ and (a) $\widetilde {Ca} = 0.01$, (b) $\widetilde {Ca} = 0.001$ at different times as indicated in the figure inset. Panel (c) depicts the growth rate of the perturbations with time.

Figure 18. For an initial perturbation that triggers the interface mode, the streamlines for the flow and the density field at time $t=135$. The parameter space is the same as in figure 16, $Pe = 5$, $\beta = 2$, $\widetilde {Ca} = 0.001$ and $k = 0.1$.

4.2.3. Pullers with a positive gradient

The instability of a suspension of pullers is special in the sense that although the interfacial and active stress strive to dampen the interface and density perturbations respectively, the resultant velocity field ends up amplifying the perturbations to the system. Figure 19(a) displays the flow streamlines which indicate that the swimmers accumulate near the minima of the interface. This is aided by the interfacial stress as seen from figure 19(d), which varies linearly with the concentration of swimmers at the interface. On the other hand, the active stress generated by the suspension of pullers amplifies the interface perturbation, as seen from figure 19(c). The instability of a suspension of pullers is thus dependent on both the interfacial and active stresses generated within the system.

Figure 19. Flow streamlines and puller density field for an imposed attractant gradient along the positive $z$-direction, for $Pe = 5$, $\beta = -20$, $k = 0.1$ and $\widetilde {Ca} = 0.001$ at time $t = 125.5$.

The effects of surface tension play a role in stabilizing the perturbations and bringing about the nonlinear saturation of the growth rate as seen from figure 20(c). At high values of $\widetilde {Ca}$, the weak effects of surface tension (figure 20a) leads to the formation of waves on the interface with the swimmers clustering at the troughs associated with the interface height. However, the time evolution is limited to numerical difficulties encountered in resolving the steep growth associated with the extremum that develops in the problem.

Figure 20. Nonlinear evolution of interface and puller density field for an imposed attractant gradient along the positive $z$-direction ($g=+1$), for $Pe = 5$, $\beta = -20$, $k = 0.1$ and (a) $\widetilde {Ca} = 0.01$, (b) $\widetilde {Ca} = 0.001$ at different times as indicated in the figure inset. Panel (c) depicts the growth rate of the perturbations with time.

5. Conclusion and discussion

In this article we have detailed a new mode of instability associated with a free interface that arises in an active suspension of swimmers. We have termed this the interface mode owing to its destabilizing effect on the system even in the absence of perturbations to the swimmer density field (as elaborated in the toy problem formulated in Appendix C). In § 3 we start by formulating a linear stability theory for the system. Through a long-wave stability calculation in § 3.1 we show an instability in suspensions of pullers, arising from the coupling between the interface and density modes. In particular, through the toy problem formulated in Appendix D we pin down the source of the instability to be the jump in normal stress at the interface created by perturbations to the swimmer density field. The stability of the system to finite wavenumber perturbations in discussed in § 3.2. The observed instabilities are found to be sensitive to the direction of the imposed attractant gradient with respect to the interface. In § 4, a nonlinear long-wave theory governing the evolution of the film is derived. The numerical solution of the nonlinear system is used to validate the growth rates estimated from the linear stability calculation. In addition, the simulations allow us to ascertain and verify the mechanism associated with the instabilities found in the linear stability calculation.

The instability of the system through the interface and density mode analysed in this paper originates from the normal stress difference that arises due to the preferential swimming of micro-swimmers along the local attractant gradient. The current work assumes a linear variation of the attractant concentration along the vertical direction. However, an accurate modelling of the system requires accounting for the role of convective transport of the attractant and the effects of secretion or consumption of attractant by micro-swimmers, as shown below in the scalar transport equation for the attractant concentration field

(5.1)$$\begin{gather} \frac{\partial c}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left[ \boldsymbol{u} c - D_a \boldsymbol{\nabla} c \right] = \alpha n. \end{gather}$$

Here, $c$ is the attractant concentration, $D_a$ depicts the diffusivity of the attractant arising from the thermal fluctuations in the system and $\alpha$ is the rate of consumption or secretion of the attractant by an individual swimmer. The approximation of the attractant concentration profile as a linear gradient is valid in the limit of the attractant diffusivity being much larger than the swimmer diffusivity $( D_a \gg D )$ and the time scale associated with the attractant secretion or consumption being much smaller than the attractant diffusive time scale $(H^2 / D_a \ll c_0 / ( \alpha \langle n \rangle )$. A recent manuscript by us (Murugan & Roy 2022) details the role of attractant transport on the creation of hydrodynamic instabilities in an active suspension. However, the addition of an attractant density field, evolving with the dynamics of the film in the current continuum framework adds a new level of complexity that is beyond the scope of the current work.

The dynamics of the attractant transport represents one possible mechanism of influencing the orientation distribution of swimmers in the suspension and consequently the active stress in the suspension. In addition, the orientation distribution is directly influenced by the viscous torques present in the system, which alters the effect of the active stress in the suspension. To fully account for the effects of shear rotation, it is necessary to solve the full kinetic equation as given in Kasyap & Koch (2014)

(5.2)\begin{equation} \mathcal{H}^{{-}1} \left\{ \frac{\partial \varOmega}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} \left[ ( \boldsymbol{p} + \boldsymbol{u} ) \varOmega \right] + \boldsymbol{\nabla}\boldsymbol{\cdot} \left( \dot{\boldsymbol{p}} \, \varOmega \right) \right\} + \left\{ \frac{\varOmega}{\tau^*} - \frac{1}{4 {\rm \pi}} \int \frac{\varOmega}{\tau^*} \, {\rm d}\boldsymbol{p} \right\} = 0, \end{equation}

where $\mathcal {H} = H / U_s \tau _0$ depicts a Péclet number based on the swimming speed of the organism. Kasyap & Koch (2014) have analytically solved the above kinetic equation through a perturbation expansion in the limit of $\mathcal {H}^{-1} \ll 1$ and obtained the following leading-order modification to the active stress:

(5.3)\begin{equation} \boldsymbol{T^{B}} = \beta \left\{ \frac{\xi}{6} \left( \boldsymbol{g g} - \frac{\boldsymbol{I}}{3} \right) n - \frac{1}{5} \mathcal{H}^{{-}1} n \left( \nabla \boldsymbol{u} + \boldsymbol{\nabla} \boldsymbol{u}^{{{\dagger}}} \right) \right\}. \end{equation}

In the above approximation of the active stress, $\beta / \mathcal {H} = {O} (\langle n \rangle L^3 ({U_s \tau _0 }/{L}))$ (Kasyap & Koch 2014). The current work is thus restricted to dilute suspensions $( \langle n \rangle L^3 \ll 1 )$ of active swimmers, wherein the effects of shear rotation may be neglected.

While the current work has addressed the evolution of the interface evolution through the active stress generated by the swimmers, the behaviour of a force dipole near a wall or interface is non-trivial, as shown by Crowdy et al. (2010), Spagnolie & Lauga (2012). In addition, in the case of $g=+1$, the clustering of swimmers near the interface may give rise to altered surface properties of surface tension and viscosity. Furthermore, the clustering of swimmers could give rise to a jammed state as seen in Dixit & Homsy (2013) and Vella, Aussillous & Mahadevan (2004), wherein the behaviour of the interface is akin to an elastic sheet. The inclusion of these effects associated with the confinement in the formulation of the active stress would lead to better predictions involving both the stability calculation and the nonlinear evolution.

The numerical solutions from the long-wave theory developed in § 4 for the evolution of the film reveals a behaviour (see figures 15c and 17b) that bears a resemblance to the rupture problem driven by Marangoni stresses as seen in Oron & Rosenau (1994), Joo, Davis & Bankoff (1991) and Kabova et al. (2006). In the limit of $Pe \ll 1$, the swimmer density field at leading order is uniform across the system domain and it can be shown that the nonlinear long-wave equations for the dimensional height evolution of the film reduce to the following equation:

(5.4)\begin{equation} h_{t} + \frac{\partial}{\partial x} \left[ \frac{h^3}{3 \mu} \frac{\partial }{\partial x} \left( \gamma h_{xx} - \frac{3\beta\mu U_0}{32}\log \frac{H}{h} \right) \right] = 0,\end{equation}

which greatly resembles the nonlinear equation derived by Burelbach, Bankoff & Davis (1988) (see also Zhang & Lister 1999) for the van der Waals rupture of a thin film

(5.5)\begin{equation} h_t + \frac{\partial}{\partial x} \left[ \frac{h^3}{3 \mu} \frac{\partial}{\partial x} \left( \gamma h_{xx} - \frac{A^*}{6 {\rm \pi}h^3} \right) \right] = 0.\end{equation}

Thus, within the ambit of our present modelling approach, active-stress-driven instabilities of interfacial flows could act as potential candidates for initiating rupture of thin films, with a swimming analogue of a disjoint pressure. We must add a caveat that any serious efforts in studying the rupture of a micro-swimmer-laden thin film must acknowledge that a continuum version of active stress will be suspect in the vicinity of a rupture event.

Setting up a system involving a linear variation of the attractant is an arduous task owing to the hydrodynamic fluctuations in the active suspension. The set-up by Cheng et al. (2007) is appropriate for studying chemotaxis in a channel flow. Experiments involving pusher- and puller-type swimmers such as E. coli and C. reinhardtii, respectively, using a set-up similar to Sokolov et al. (2009), probably involving phototaxis for suspensions of pullers, would provide further insight into the instability caused through the coupling between the active and interfacial stresses in the film.