Localization and bistability of bioconvection in a doubly periodic domain

We consider a two-dimensional model convection system where the velocity field $\bm{u}(x,y,t)=(u_{x},u_{y})$ is coupled with the density deviation $m(x,y,t)$ from an equilibrium density profile $m_{0}(y)=\sin(y)$ in a doubly periodic domain $(x,y)\in\Omega\equiv[0,L_{x}]\times[0,2\pi]$. We assume that $\bm{u}$ is incompressible ($\nabla\cdot\bm{u}=0$) and thus the velocity field $\bm{u}$ is represented by a stream function $\Psi(x,y,t)$: $(u_{x},u_{y})=(\partial_{y}\Psi,-\partial_{x}\Psi)$.

As an extended model of fluid convection, we consider the following nondimensional equations for $m$ and the vorticity $\omega(x,y,t)=\partial_{x}u_{y}-\partial_{y}u_{x}$:

	$\displaystyle\frac{\partial\omega}{\partial t}-J(\Psi,\omega)$	$\displaystyle=Pr(\Delta\omega-Ra\partial_{x}m),$		(1)
	$\displaystyle\frac{\partial m}{\partial t}-J(\Psi,m)$	$\displaystyle=\Delta m-Pe\partial_{y}m+\partial_{x}\Psi\partial_{y}m_{0},$		(2)

where $J(\Psi,f)$ ($f$ is $\omega$ or $m$) is defined as

$$J(\Psi,f)\equiv(\partial_{x}\Psi)(\partial_{y}f)-(\partial_{x}f)(\partial_{y}\Psi).$$

(3)

Derivations of the model with its physical background are provided in the next subsection [Sec.II.2]. Here, three nondimensional parameters, $Ra,Pe,Pr$ represent the Rayleigh number, Prandtl number $Pr$, and Péclet number, respectively. The Péclet number represents nondimensional particle self-propulsive velocity, and when the Péclet number approaches zero, Eqs.(1)-(2) recover the Rayleigh-Bénard thermal convection, in which the last term in Eq.(2) is simply reduced to $\partial_{x}\Psi$, because of the linear profile of the equilibrium density. The stream function $\Psi$ satisfies the Poisson equation on the doubly periodic domain $(x,y)\in\Omega$ as follows:

$\displaystyle-\Delta\Psi$

$\displaystyle=\omega.$

(4)

We readily find that Eqs.(1)-(2) have a quiescent state as a trivial solution, $\omega=m=0$, which is a steady solution for any $(Ra,Pr,Pe)$.

In this subsection, we derive our physical model, Eqs. (1)-(2), and discuss the underlying physical backgrounds. We follow a standard bioconvection model for a bulk motion[2, 24, 3], but with a doubly periodic boundary condition. By considering an equilibrium density profile as an independent model parameter, we naturally extend the usual bioconvection model to a general boundary condition.

Following well-known continuum description for a dilute solution of self-propelled particles, we introduce the density field $\tilde{m}(x,y,t)$ and the velocity field for the dilute solution $\bm{u}(x,y,t)=(u_{x}(x,y,t),u_{y}(x,y,t))$, where we assume the incomprehensibility condition, $\partial_{x}u_{x}+\partial_{y}u_{y}=0$. As the simplest case of bioconvection, we follow Fick’s law with a constant diffusion for the density flux, $\bm{J}$, as

$\displaystyle J(\tilde{m})$

$\displaystyle=\tilde{m}\bm{u}+\tilde{m}\bm{V}-D\bm{\nabla}\tilde{m}.$

(5)

Here, we introduce upward locomotion of the self-propelled motion with $\bm{V}=V\bm{\hat{y}}$, where $\bm{\hat{y}}$ is the unit vector in the $y$ direction and $V$ and $D$ are constant values. The upward motion of self-propelled particles is considered to be essential for bioconvective pattern, and in experiments these biased motions usually originate from the cellular tactic behaviors, e.g. phototaxis, gyrotaxis, and chemotaxis.

Then, we define an equilibrium density profile $m_{0}(y)$ as the steady solution with $\bm{u}=0$. We then introduce the density deviation from the equilibrium profile as $m=\tilde{m}-m_{0}$. The density flux for $m$ is then obtained by $\delta J\equiv J(\tilde{m})-J(m_{0})$, and we have

$\displaystyle\delta J$

$\displaystyle=m\bm{u}+mV\bm{\hat{y}}+m_{0}\partial_{y}u_{y}-D\bm{\nabla}m.$

(6)

Hence, the time evolution of $m$ follows the continuum equation, $\partial_{t}m=-\bm{\nabla}\cdot\delta J$, which is explicitly given by

$\displaystyle\frac{\partial m}{\partial t}+\bm{u}\cdot\bm{\nabla}m$

$\displaystyle=-V\partial_{y}m-u_{y}\partial_{y}m_{0}+D\Delta m,$

(7)

where $\Delta=\partial_{x}^{2}+\partial_{y}^{2}$ is Laplacian in two dimensions.

For the fluid motion, we consider the two-dimensional Navier-Stokes equation with Boussinesq approximation. The velocity $\bm{u}(x,y,t)$ and the pressure $p(x,y,t)$ then satisfy

$\displaystyle\rho_{s}\left(\frac{\partial\bm{u}}{\partial t}+\bm{u}\cdot\bm{\nabla}\bm{u}\right)$

$\displaystyle=-\bm{\nabla}p+\mu\Delta\bm{u}-mv(\rho_{p}-\rho_{s})g\bm{\hat{y}}.$

(8)

Here, $\mu$ is the dynamic viscosity and assumed to be a constant. $\rho_{s}$ and $\rho_{p}$ are the density of solvent and particles, respectively, and in the context of bioconvection, we usually consider $\rho_{p}>\rho_{s}$. The gravity constant $g>0$ and the particle volume $v$ are also assumed to be constant.

We finally derive our equation of motions Eqs.(1)-(2) by introducing the stream function $\Psi$. Three nondimensional parameters, Rayleigh number $Ra$, Peclet number $Pe$ and Prandtl number $Pr$, are then introduced with physical parameters as follows:

	$\displaystyle Ra$	$\displaystyle=\frac{||m_{0}||_{\infty}v(\rho_{p}-\rho_{s})gL_{y}^{3}}{8\pi^{3}D\mu},$		(9)
	$\displaystyle Pe$	$\displaystyle=\frac{L_{y}V}{2\pi D},$		(10)
	$\displaystyle Pr$	$\displaystyle=\frac{\mu}{\rho D},$		(11)

where $L_{y}$ is the length of domain in $y$. We have assumed the length scale $L=L_{y}/2\pi$, the time scale $T=L^{2}/D$, and the density scale as $||m_{0}||_{\infty}$. Our model includes Rayleigh–Bénard convection, which is the most fundamental model for thermal convection, as the case for $Pe=0$ and $dm_{0}/dy=$constant. If a non-slip or free-slip boundary condition is imposed on the bottom and top boundaries of the fluid region, the equilibrium density profile is then accordingly determined. In this sense, our extended model captures the bulk behaviour of bioconvection.

Focusing on horizontal accumulation of particles, we introduce vertically averaged density $n(x,t)\equiv\langle m\rangle$, where the bracket $\langle\rangle$ is average in $y$, defined as

$\displaystyle\langle\alpha\rangle\equiv\frac{1}{2\pi}\int_{0}^{2\pi}\alpha(x,y,t)dy.$

(12)

for a field variable $\alpha(x,y,t)$. By integrating Eq.(2) in the $y$ direction, we find

$\displaystyle\frac{\partial n}{\partial t}$

$\displaystyle=\partial_{x}^{2}n+\partial_{x}\langle(m_{0}+m)u_{x}\rangle.$

(13)

This equation suggests that the particle accumulation can occur only when the second term of Eq.(13) has a nonzero contribution.

To physically interpret this nonlinear effect, we introduce a horizontal net density flux as $F_{x}\equiv-\langle(m_{0}+m)u_{x}\rangle$. As shown in Fig.1(a), let us consider the situation where the particle concentration is denser in the upper region [also illustrated by the green shades in Fig. 1(b)]. To generate the localized convection pattern from the horizontally homogeneous trivial state, nonzero horizontal net flux must be produced as indicated by magenta arrows in Fig.1(b). Since the sign of $F_{x}$ is dominated by the sign of $u_{x}$ in the region of higher $m+m_{0}$ values, the horizontal fluid velocity $u_{x}$ needs to be generated in the same direction as the net density flux at the upper region, leading to an incoming flow at the point of the particle condensation and the associated outgoing flow at the lower region due to the fluid incompressibility. In the end, these mechanical couplings result in the downward jet accompanied by the particle condensation with a circulatory fluid motion as in the blue arrows in Fig.1(b).

These physical mechanisms imply that the trivial solution could be destabilized by the nonlinear coupling and then the system forms a large-scale convective pattern. We therefore focus on the large-scale coupling to capture the dominant nonlinear effects of the system. In doing so, we introduce a minimal system by truncating the higher Fourier modes in the vertical direction. We retain only the lowest three Fourier modes which are proportional to $\exp(ily)$ with $l=0,\pm 1$ as

	$\displaystyle\omega(x,y,t)$	$\displaystyle=A_{-1}(x,t)\exp(-iy)+A_{0}(x,t)+A_{1}(x,t)\exp(iy),$		(14)
	$\displaystyle m(x,y,t)$	$\displaystyle=B_{-1}(x,t)\exp(-iy)+B_{0}(x,t)+B_{1}(x,t)\exp(iy),$		(15)
	$\displaystyle\Psi(x,y,t)$	$\displaystyle=C_{-1}(x,t)\exp(-iy)+C_{0}(x,t)+C_{1}(x,t)\exp(iy),$		(16)

where $i=\sqrt{-1}$ is the imaginary unit. To validate our truncation model, we have also numerically confirmed that the qualitative behaviors of the system are unchanged even in the presence of the higher modes, $l=0,\pm 1,\cdots,\pm 64$. Henceforth, in this paper, we investigate the minimal large-scale model with Eqs. (14)-(16).

Since $\omega$ and $m$ are variables with real value, $A_{0}$ and $B_{0}$ must also be real. Similarly, $A_{1}^{\ast}=A_{-1}$ and $B_{1}^{\ast}=B_{-1}$ hold, where the asterisk indicates complex conjugate. From the Poisson equation (4), the stream function $\Psi$ satisfies

$\displaystyle-(\partial_{x}^{2}-l^{2})C_{l}$

$\displaystyle=A_{l},$

(17)

for $l=0$, and $1$.

We first expand $A_{l}$, $B_{l}$, and $C_{l}$ ($l=0,\pm 1$) by horizontal Fourier series as

	$\displaystyle A_{l}(x,t)$	$\displaystyle=\sum_{k=-M}^{M}\tilde{A}_{k,l}(t)\exp(ikx),$		(18)
	$\displaystyle B_{l}(x,t)$	$\displaystyle=\sum_{k=-M}^{M}\tilde{B}_{k,l}(t)\exp(ikx),$		(19)
	$\displaystyle C_{l}(x,t)$	$\displaystyle=\sum_{k=-M}^{M}\tilde{C}_{k,l}(t)\exp(ikx),$		(20)

where we choose a truncation wavenumber as $M=64$ to sufficiently resolve horizontal structures. Hence, our model equation, Eqs. (14) - (16), are discretized into the following equations:

	$\displaystyle\left(\frac{d}{dt}+Prk^{2}\right)\tilde{A}_{k,0}$	$\displaystyle=-ikPrRa\tilde{B}_{k,0}+\sum_{p,q}qk\frac{\tilde{A}_{p,q}\tilde{A}_{k-p,-q}}{p^{2}+q^{2}},$		(21)
	$\displaystyle\left(\frac{d}{dt}+Pr(k^{2}+1)\right)\tilde{A}_{k,1}$	$\displaystyle=-ikPrRa\tilde{B}_{k,1}+\sum_{p,q}(qk-p)\frac{\tilde{A}_{p,q}\tilde{A}_{k-p,1-q}}{p^{2}+q^{2}},$		(22)
	$\displaystyle\left(\frac{d}{dt}+k^{2}\right)\tilde{B}_{k,0}$	$\displaystyle=ik\frac{\tilde{A}_{k,1}+\tilde{A}_{k,-1}}{2(k^{2}+1)}+\sum_{p,q}qk\frac{\tilde{A}_{p,q}\tilde{B}_{k-p,-q}}{p^{2}+q^{2}},$		(23)
	$\displaystyle\left(\frac{d}{dt}+(k^{2}+1)+iPe\right)\tilde{B}_{k,1}$	$\displaystyle=ik\frac{\tilde{A}_{k,0}}{k^{2}}+\sum_{p,q}(qk-p)\frac{\tilde{A}_{p,q}\tilde{B}_{k-p,1-q}}{p^{2}+q^{2}}.$		(24)

Here, we used the relation, $C_{k,l}=(k^{2}+l^{2})^{-1}A_{k,l}$, which follows from Eq.(17).

Notably, the system has a reflection symmetry in $x$, as Eqs.(1) and (2) are unchanged by a reflection $R_{x}:(\omega(x,y,t),m(x,y,t))\rightarrow(-\omega(-x,y,t),m(x,y,t))$, which yields the relation between the horizontal Fourier coefficients $\tilde{A}_{k,l}=-\tilde{A}_{-k,l}$ and $\tilde{B}_{k,l}=\tilde{B}_{-k,l}$. Since the real condition of $\omega$ and $m$ read $\tilde{A}_{k,l}=\tilde{A}_{-k,-l}^{\ast}$ and $\tilde{B}_{k,l}=\tilde{B}_{-k,-l}^{\ast}$, we obtain the relation, $\tilde{A}_{k,-l}=-\tilde{A}^{\ast}_{k,l}$ and $\tilde{B}_{k,-l}=\tilde{B}_{k,l}^{\ast}$, which, for example, guarantees that the first term in the right-hand side of Eq.(23) is real.

By this symmetry, $\tilde{A}_{k,0}$ and $\tilde{B}_{k,0}$ are pure imaginary and real numbers, respectively, and the time evolution of our system, Eqs. (21)-(24), is determined by a set of Fourier modes with non-negative $k$ indices, leading to a dynamical system of a state $\bm{X}(t)$ with $(M+1)\times 6$ real variables.

In our numerical computations, we have utilized the pseudospectral method for computing nonlinear terms complemented by a Fast Fourier Transform (FFT). In time stepping, the linear terms have been dealt as the integral factor and we have solved the nonlinear terms by the 4th-order Runge-Kutta method.

To compute stationary solutions, we have applied the GMRES-Newton method, in which the linear space is approximated in a Krylov subspace. Iterations of Newton’s method have been terminated when the time derivatives of the Fourier modes are sufficiently small. More precisely, we employed the condition, $||d\bm{X}/dt||_{2}<10^{-5}$, where $||\bullet||_{2}$ indicates the Euclidean norm. The stationary solutions in different parameters are continuously obtained by the arclength continuation method.

We then proceed to analyze the nonlinear regime, focusing on smaller Péclet numbers where the trivial solution becomes linearly unstable. To evaluate the impact of $Pe$ on the convection structure and its stability, we performed a bifurcation analysis with the horizontal scale of the fluid region fixed as $L_{x}=8\pi$, following the existing literature on the thermal convection and forced Navier-Stokes flow problems in a doubly periodic domain [24, 25, 9, 10]. The critical Rayleigh number at $Pe=0$ is then obtained from Eq. (29) as $Ra=\sqrt{2k_{s}^{2}(1+k_{s}^{2})^{3}}\approx 0.388$. Since our current focus is on the nonlinear regime, we choose the Rayleigh number to be larger than this critical value and set $Ra=0.5$. We also set the Prandtl number $Pr=2.5$ as a reasonable value in the experimental setup of bioconvection of Euglena [26]. This choice of value is also reasonable in study of Paramecium by Taheri & Bilgen [24] in which the nondimensional parameter ranges are estimated as $Pr=0.22-2$ and $Pe=0.07-1.54$.

We first show a bifurcation curve in Fig.2(a), where the max norm of $n$, $||n||_{\infty}$, is plotted against the Péclet number $Pe$. The stable and unstable steady solutions are colored red and blue, respectively. At $Pe=0$, a nonlinear stable solution exists in addition to the unstable trivial solution, and a unstable nonlinear solution bifurcates from the trivial solution at the critical Péclet number, $Pe=Pe_{c}$. The stable and unstable nonlinear solutions collapse at $Pe_{SN}\approx 4.27$, where a saddle-node bifurcation occurs. In the region of $Pe\in(Pe^{c},Pe_{SN})\approx(0.868,4.27)$, there exist two nonlinear solutions, and we hereafter call the stable solution as UB (upper branch) solution denoted by the state vector $\bm{X}_{\rm UB}$ and the unstable solution as LB (lower branch) solution, which we denoted by $\bm{X}_{\rm LB}$.

In Fig.2(b,c), we show density deviation fields, $m(x,y)$, for the upper and lower branches at $Pe=3.8$, which correspond to the red and blue circles in Fig.2(a), respectively. In both cases, the density deviation field is spatially localized. Note that the system possesses a horizontal translational symmetry, and we choose a solution that is localized at the middle of the domain. The differences between the two steady solutions are discussed detailed in Sec. IV.2.

Beyond the saddle-node bifurcation $(Pe>Pe_{SN})$, we numerically confirmed that there is only a stable trivial solution. Since the Péclet number represents the nondimensional self-propulsive velocity, these results suggest that the self-motility could induce the bistable structure in the convection system.

We then focus on the spatial structure of the fluid and particle density fields. To characterize steady solutions, we consider a vertical average of the upward velocity and a similar vertical average of the density deviation, which are denoted by $U_{y}\equiv\langle u_{y}(x,y)\rangle=-\partial_{x}C_{0}$ and $n(x)=\langle m(x,y)\rangle$, respectively. We plot these values for the nonlinear stable solution at $Pe=0$ in Fig.3(a) and Fig.3(b), as well as the two complex modes of the vorticity and density deviation fields [see Eqs. (14)-(15)], $A_{1}(x)$ and $B_{1}(x)$, in Fig.3(c) and Fig.3(d). The real and imaginary parts are shown as solid and broken lines, respectively.

From these plots, we find that the field structure is horizontally separated into two region; one is characterized by the upward fluid velocity and the dilute density field, and the other has the opposite properties. This anticorrelation between the fluid and density fields agrees with the schematics of the nonlinear coupling in Fig.1. Also, this result is consistent with experimental observations of bioconvection [27, 5, 28, 3]. We observe that the spatial separation is visible more clearly for the density field, as seen in the kink-like structures of the figure. The complex amplitudes in Fig.3(c,d) exhibit peaky structures around the kinks seen in the plots of $U_{y}$ and $n$.

We then examine the same stable nonlinear solution but with an increased Péclet number. In Fig.4, we show the similar plots to characterize the spatial structure at $Pe=3.8$.

As in Fig. 3, we find the large-scale spatial separation with the anticorrelation between the velocity and density fields. Compared with the density concentration at $Pe=0$, the density field $n(x)$ in Fig.4(b) exhibits a further localized structure. To quantify the width between the kink and anti-kink, we introduce a size of the concentrated region as $L_{n}=\frac{\min n}{\min n-\max n}L_{x}$. By numerically evaluating the values of $L_{n}$ in different $Pe$, we found that as $Pe$ increases $L_{n}$ monotonically decreases from $L_{n}=4\pi=0.5L_{x}$ at $Pe=0$ to $L_{n}\approx 0.204L_{x}$ at $Pe=Pe_{\rm SN}$.

We then proceed to study the LB solution, which forms the unstable lower branch in Fig.2. We examine the spatial structure of $X_{\rm LB}$ at $Pe=3.8$ and the results are summarized in Fig.5. As shown in Fig.5(a,b), the LB solution exhibits spatially-localized profiles both in the vertical velocity and density concentration than the UB solution, although the amplitudes of the profile are slightly suppressed.

We also show the complex amplitudes in the vorticity and density deviation fields in Fig.5(c,d). Notably, the $B_{1}$ amplitude does not possess the peaks around the kinks of the localized field [Fig.5(d)]. This difference are clearly shown in the plot of $|B_{1}|$ in Fig.6(a), as seen in the double-peak structure for the UB solution (red) and single-peak structure for the LB solution (blue).

To further quantify these differences, we introduce the distance between the two peaks in $|B_{1}|$ as $L_{B_{1}}$. Hence, $L_{B_{1}}=0$ when two peaks merge, and we found $L_{B_{1}}=0$ for all LB solutions. The value of $L_{B_{1}}$ is then calculated for the UB solution by estimating the peak position using the three-point Lagrange interpolation method, the points of the grid being equally separated into 256 parts. As shown in Fig. 6, the value of $L_{B_{1}}$ for $\bm{X}_{\rm UB}$ decreases monotonically to reach zero around $Pe\approx Pe_{\rm SN}$.

In the parameter regime discussed earlier in this section, we have observed the bistable nature of the system. We then further examine the global stricture of the dynamical system.

By analyzing the time evolution around the unstable UB solution, we numerically observed some orbits stay longer at the neighborhood. We then chased the generated orbits and using the bisection method (details are explained later), we have numerically detected the basin boundary by which the final asymptotic state is separated. In Fig.7(a), we show schematic of the basin boundary and the global structure of the steady solutions. While the orbits in the state space are attracted to the trivial solution or the UB solution, the LB solution behaves as an edge state [20, 23], which separates the attracting state of orbits. As schematically shown in the figure, the unstable manifold ($W_{u}$) connects to each stable fixed point, and the stable manifold ($W_{s}$) forms a basin boundary between the trivial solution and the UB solution.

Indeed, we have numerically obtained these structures and the results are shown in Fig.7(b), where the computed trajectories are projected onto a two-dimensional state space that is spanned by the norms$(||X||_{2},||X-X_{\rm LB}||_{2})$. The two stable states are shown in the colored circles: the trivial solution marked by a green circle at the upper-left corner and the UB solution marked by a red circle at the right-upper corner. The two representative trajectories are shown in different colors, and the green and red curves show orbits reaching the trivial and UB solutions, respectively, corresponding the schematic trajectories in Fig.7(a). The dotted lines represent the basin boundary that is numerically obtained by the bisection method.

We now briefly describe the bisection method that is used to detect the basin boundary. We first consider a state as a linear interpolation of the two steady solution as $\bm{Y}_{\alpha}=\alpha\bm{X}_{\rm UB}(\alpha\in[0,1])$. Our aim is then to find $\alpha^{c}\in[0,1]$, around which the eventual states are separated [Fig.7(a)]. The numerical results show that the orbits stay at the neighborhood of the LB solution, $\bm{X}_{LB}$, and depart exponentially in time, as shown in the schematic orbits [Fig.8(a)] and in the projection [Fig.8(b)]. Further details are shown in Fig.8, where we show several orbits near the basin boundary. In Fig.8(a) we present a set of orbits that finally reach the trivial solution (green) and the UB solution (red), and the basin boundary is then obtained as the broken line. The plots in Fig.8(b) show the relative distance to the LB solution for the same trajectories and we find that the states initially located nearer to the basin boundary stay longer at the neighborhood of the LB solution.