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 Reference Subramanian and Koch2009)

(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 (Reference Subramanian and Koch2009)

(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 (Reference Subramanian, Koch and Fitzgibbon2011) 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. (Reference Chen, Ford and Cummings2003) and Rivero et al. (Reference Rivero, Tranquillo, Buettner and Lauffenburger1989) 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 Reference Kasyap and Koch2014). 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 Reference Batchelor1970)

(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 (Reference Kasyap and Koch2014). 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 Reference Kasyap and Koch2012). We would like to mention the presence of a typographical error in the definition of $\beta$ in Murugan & Roy (Reference Murugan and Roy2022). 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. Reference Guasto, Johnson and Gollub2010). 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. Reference Horvat, Pannuri, Romeo, Dogsa and Stopar2019). 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 Reference Brenner, Levitov and Budrene1998). 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. Reference Sokolov, Goldstein, Feldchtein and Aranson2009) 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. Reference Guasto, Johnson and Gollub2010) 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 Reference Subramanian and Koch2009; Guasto et al. Reference Guasto, Johnson and Gollub2010). 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. Reference Sokolov, Goldstein, Feldchtein and Aranson2009), 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. (Reference Sokolov, Goldstein, Feldchtein and Aranson2009) 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 (Reference Kasyap and Koch2012) 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 Reference Leal2007). One must note that suspensions of passive non-Brownian particles also exhibit normal stresses (Guazzelli & Pouliquen Reference Guazzelli and Pouliquen2018). 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 Reference Dhas and Roy2022). 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 (Reference Kasyap and Koch2012). In Appendix B we detail the effects of a wall with slip on the stability calculation by Kasyap & Koch (Reference Kasyap and Koch2012), 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 (Reference Smith1990). 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 (Reference Kasyap and Koch2012). 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 (Reference Kasyap and Koch2012) 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 (Reference Kasyap and Koch2012, Reference Kasyap and Koch2014). 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 (Reference Kasyap and Koch2014) 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 (Reference Trefethen2000) 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 (Reference Kasyap and Koch2012) 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 (Reference Oron and Rosenau1994). 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 (Reference Kasyap and Koch2012), 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.