Effect of rotation on anisotropic scattering suspension of phototactic algae

Bioconvection is a widespread phenomenon in fluid dynamics that involves the formation of patterns in suspensions of living microorganisms, such as algae and bacteria Platt (1961); Pedley and Kessler (1992); Hill and Pedley (2005); Bees (2020); Javadi et al. (2020). The term "bioconvection" was first introduced by Platt in 1961. Typically, these microorganisms have a higher density than the surrounding fluid on a small scale and collectively swim in an upward direction. However, there are cases in natural environments where density differences are not necessary for pattern formation, and the microorganisms exhibit different behaviours. Non-living microorganisms do not exhibit pattern formation behaviour. Microorganisms respond to external stimuli and exhibit behavioural changes in their swimming direction, known as taxis. Examples of taxis include gravitaxis, which is influenced by gravitational acceleration, gyrotaxis, which is influenced by both gravitational acceleration and viscous torque when microorganisms have bottom-heavy structures, and phototaxis, which is the movement of microorganisms towards or away from a light source. This article focuses specifically on the effect of phototaxis.

The phenomenon of bioconvection patterns in algal suspensions has been the subject of extensive research, highlighting the significant role of light in shaping fluid dynamics and cellular distribution Wager (1911); Kitsunezaki, Komori, and Harumoto (2007). When exposed to intense light in well-mixed cultures of photosynthetic microorganisms, dynamic patterns can either persist or undergo disruption. Strong light can hinder the formation of stable patterns, and it can also disrupt pre-existing  Kessler (1985); Williams and Bees (2011); Kessler (1989) The response of algae to light of varying intensities is a key determinant in the modulation of bioconvection patterns. Moreover, the interaction of light with microorganisms involves processes of light absorption and scattering. Algal light scattering can be characterized as either isotropic or anisotropic, the latter further classified into forward and backward scattering based on cell properties such as size, shape, and refractive index. Algae, due to their size, predominantly exhibit forward scattering of light within the visible wavelength range.

The phototaxis model proposed by Rajput is being employed to investigate the bioconvection system. In this study, the suspension of phototactic microorganisms is assumed to undergo rotation around the z-axis at a constant angular velocity, while being illuminated from above by collimated flux. Given that many motile algae rely on photosynthesis for their energy needs, they exhibit distinct phototactic behavior. To ensure a comprehensive understanding of their response to light, it is essential to examine the phototaxis model in the context of a rotating medium. Unlike previous models that considered scattering effects in an isotropic manner, our approach incorporates anisotropic scattering. This allows for a more precise representation of the swimming behaviour of microorganisms as they react to light stimuli. By considering the influence of anisotropic scattering, we aim to enhance the accuracy of the model and provide valuable insights into the intricate dynamics of the bioconvection phenomenon.

In this experimental setup, the mathematical formulation of the problem involves considering the behavior of a forward scattering algal suspension confined between two parallel boundaries in the (y, z) plane. The suspension is assumed to be uniformly illuminated from above by collimated flux. The goal is to study the behavior of the suspension under the influence of light.

To analyze the system, a realistic phototaxis model is employed, which takes into account the forward scattering of light by the algae. The equilibrium solution of the system is determined by considering a steady-state where the flow velocity is zero, indicating no fluid motion. In this state, the balance between phototaxis and diffusion leads to the formation of a horizontal sublayer where cells accumulate, dividing the suspension into two distinct regions.

The critical intensity $G_{c}$ plays a crucial role in determining the position of the sublayer. Above the sublayer ($G>G_{c}$), the light intensity is sufficient to suppress cell movement, while below the sublayer ($G<G_{c}$), cells exhibit upward motion towards the sublayer due to phototaxis.

The lower region beneath the sublayer is considered unstable and has the potential to exhibit fluid motion and penetrate the stable upper region through penetrative convection if the system becomes unstable. Penetrative convection is a phenomenon observed in various convection problems.

Over the years, researchers have made significant advancements in modeling phototaxis and bioconvection, employing various approaches such as linear stability theory and numerical simulations. For instance, Vincent and Hill Vincent and Hill (1996) investigated the initiation of phototactic bioconvection and identified non-stationary and non-oscillatory disturbance modes. Ghorai and Hill Ghorai and Hill (2005) conducted a two-dimensional simulation of bioconvective flow patterns based on Vincent and Hill’s model, albeit without considering algae scattering. Ghorai et al. Ghorai, Panda, and Hill (2010) examined the onset of bioconvection in an isotropic scattering suspension and observed a steady-state profile with a sublayer located at different depths due to isotropic scattering. Panda and Ghorai Panda and Ghorai (2013) conducted a linear stability analysis in a forward scattering algal suspension, observing a transition from non-oscillatory to non-stationary modes during bioconvective instability caused by forward scattering, while Panda and Singh Panda and Singh (2016) numerically explored the linear stability of PBC using Vincent and Hill’s continuum model in a two-dimensional setting. Previous studies did not account for the effects of diffuse irradiation. Panda et al. Panda et al. (2016) proposed a model that examined how diffuse radiation influenced isotropic scattering algal suspensions, revealing that diffuse irradiation significantly stabilizes such suspensions. Panda Panda (2020) further investigated the impact of forward anisotropic scattering on the onset of PBC, considering both diffuse and collimated irradiation. However, these studies did not account for the influence of oblique collimated flux. Panda et al. Panda, Sharma, and Kumar (2022) introduced oblique collimated flux and discovered that bioconvective solutions switch from non-oscillatory to overstable states and vice versa during bioconvective instability induced by oblique collimated flux. Kumar Kumar (2022) examined the effect of oblique collimated irradiation on isotropic scattering algal suspensions, revealing the existence of both stationary and overstable solutions within certain parameter ranges. Kumar Kumar (2023a) investigated the effect of rotation on the phototactic bioconvection in a non-scattering medium and found a significant stabilizing impact due to rotation. More recently, Panda and Rajput Panda and Rajput (2023) investigate the combined impact of oblique and diffuse irradiation on the onset of bioconvection.

However, there is still a need to investigate phototactic bioconvection in a rotating forward scattering algal suspension using a realistic phototaxis model. Therefore, the objective of this study is to examine bioconvective instability in such a system and provide insights into the complex dynamics of bioconvection in a light-induced pattern formation in a dilute algal suspension.

The article follows a structured approach, starting with the mathematical formulation of the problem. It derives the equilibrium solution and perturbs the governing system of bioconvection with small disturbances. The linear stability problem is formulated and solved using numerical methods. The obtained results are then presented and discussed. The article concludes by highlighting the novelty and significance of the proposed model in advancing the understanding of bioconvective instability in rotating forward scattering algal suspensions.

The research focuses on studying the behaviour of a forward scattering algal suspension within the (y, z) plane, which is confined between two parallel boundaries. The boundaries are assumed to be infinite in extent, and there is no reflection of light from the top and bottom boundaries. This setup allows for the investigation of the suspension’s behaviour under the influence of light.

The suspension is uniformly illuminated from above by collimated flux. Collimated flux refers to light that is travelling in parallel rays, indicating that the incident light rays are not scattered before reaching the suspension. This uniform illumination from above serves as the driving force for the system and influences the phototactic behaviour of the algae.

The equilibrium state of the system is characterized by the following conditions:

$$\bm{U}=0,~{}~{}~{}n=n_{s}(z),~{}~{}~{}\zeta_{s}=\nabla\times U=0,~{}~{}~{}\text{and}~{}~{}~{}L=L_{s}(z,\theta).$$

At equilibrium, the total intensity $G_{s}$ and radiative flux $\bm{q}_{s}$ can be expressed as:

$$G_{s}=\int_{0}^{4\pi}L_{s}(z,\theta)d\Omega,\quad\text{and}\quad\bm{q}_{s}=\int_{0}^{4\pi}L_{s}(z,\theta)\bm{s}d\Omega.$$

The $\bm{q}_{s}$ vector has zero components in $x$ and $y$ directions due to the independence of $L_{s}^{d}(z,\theta)$ on $\phi$. Therefore, $\bm{q}_{s}$ can be written as $\bm{q}_{s}=-q_{s}\hat{\bm{z}}$, where $q_{s}$ represents the magnitude of $\bm{q}_{s}$.

The equation governing $L_{s}$ at equilibrium is given by:

$$\frac{dL_{s}}{dz}+\frac{\kappa n_{s}L_{s}}{\nu}=\frac{\omega\kappa n_{s}}{4\pi\nu}(G_{s}(z)-Aq_{s}\nu).$$

(24)

The equilibrium intensity $L_{s}$ can be decomposed into two components: the collimated part $L_{s}^{c}$ and the diffuse part caused by scattering $L_{s}^{d}$. The equation governing the collimated component $L_{s}^{c}$ is:

$$\frac{dL_{s}^{c}}{dz}+\frac{\kappa n_{s}L_{s}^{c}}{\nu}=0,$$

(25)

with the boundary condition:

$$L_{s}^{c}(1,\theta)=L_{t}\delta(\bm{r}-\bm{r}_{0}),\quad\theta\in[\pi/2,\pi],$$

(26)

where $L_{t}$ represents the total intensity at the point source location $\bm{r}_{0}$.

Solving the governing equation for $L_{s}^{c}$ with the given boundary condition yields the expression for $L_{s}^{c}$:

$$L_{s}^{c}=L_{t}\exp\left(\int_{z}^{1}\frac{\kappa n_{s}(z^{\prime})}{\nu}dz^{\prime}\right)\delta(\bm{r}-\bm{r}_{0}).$$

(27)

The diffuse part $L_{s}^{d}$ is governed by the equation:

$$\frac{dL_{s}^{d}}{dz}+\frac{\kappa n_{s}L_{s}^{d}}{\nu}=\frac{\omega\kappa n_{s}}{4\pi\nu}(G_{s}(z)-Aq_{s}\nu),$$

(28)

with the boundary conditions:

$$L_{s}^{d}(1,\theta)=0,\quad\theta\in[\pi/2,\pi],$$

(29)

$$L_{s}^{d}(0,\theta)=0,\quad\theta\in[0,\pi/2].$$

(30)

In the basic state, the total intensity $G_{s}$ is given by:

$$G_{s}=G_{s}^{c}+G_{s}^{d}=\int_{0}^{4\pi}[L_{s}^{c}(z,\theta)+L_{s}^{d}(z,\theta)]d\Omega=L_{t}\exp\left(-\int_{z}^{1}\kappa n_{s}(z^{\prime})dz^{\prime}\right)+\int_{0}^{\pi}L_{s}^{d}(z,\theta)d\Omega,$$

(31)

where $G_{s}^{c}$ and $G_{s}^{d}$ represent the collimated and diffuse components of $G_{s}$, respectively.

Similarly, the radiative heat flux $\bm{q}_{s}$ in the basic state can be written as:

$$\bm{q}_{s}=\bm{q}_{s}^{c}+\bm{q}_{s}^{d}=\int_{0}^{4\pi}\left(L_{s}^{c}(z,\theta)+L_{s}^{d}(z,\theta)\right)\bm{r}d\Omega=L_{t}\exp\left(\int_{z}^{1}-\kappa n_{s}(z^{\prime})dz^{\prime}\right)\hat{\bm{z}}+\int_{0}^{4\pi}L_{s}^{d}(z,\theta)\bm{r}d\Omega.$$

(32)

To solve the problem, a new variable $\tau$ is defined as:

$$\tau=\int_{z}^{1}\kappa n_{s}(z^{\prime})dz^{\prime}.$$

These equations provide a set of coupled Fredholm integral equations of the second kind. Specifically, the equation for $G_{s}$ can be written as:

$$G_{s}(\tau)=e^{-\tau}+\frac{\omega}{2}\int_{0}^{\kappa}G_{s}(\tau^{\prime})E_{1}(|\tau-\tau^{\prime}|)d\tau^{\prime}+Asgn(\tau-\tau^{\prime})q_{s}(\tau^{\prime})E_{2}(|\tau-\tau^{\prime}|),$$

(33)

where $E_{1}(x)$ and $E_{2}(x)$ represent the exponential integrals of order 1 and 2, respectively.

$$q_{s}(\tau)=e^{-\tau}+\frac{\omega}{2}\int_{0}^{\kappa}Aq_{s}(\tau^{\prime})E_{3}(|\tau-\tau^{\prime}|)d\tau^{\prime}+sgn(\tau-\tau^{\prime})G_{s}(\tau^{\prime})E_{2}(|\tau-\tau^{\prime}|).$$

(34)

The coupled fractional integral equations (FIEs) can be solved using the method of subtraction of singularity, where $E_{n}(x)$ represents the exponential integral of order $n$ and $sgn(x)$ denotes the signum function.

The mean swimming direction in the basic state is given by

$$<\bm{p_{s}}>=-M_{s}\frac{\bm{q_{s}}}{q_{s}}=M_{s}\hat{\bm{z}},$$

where $M_{s}=M(G_{s})$. In the basic state, the cell concentration $n_{s}(z)$ satisfies the following equation

$$\frac{dn_{s}}{dz}-V_{c}M_{s}n_{s}=0,$$

(35)

which is accompanied by the equation

$$\int_{0}^{1}n_{s}(z)dz=1.$$

(36)

Equations (33) to (36) constitute a boundary value problem that can be solved using numerical techniques such as the shooting method.

To demonstrate the equilibrium state, we fix the values of $S_{c}=20$, $G_{c}=1$, $\kappa=0.5$, and $\omega$ at specific values of 0.1, 0.42, and 0.47. We then investigate the influence of the forward scattering coefficient $A$ by varying it within the range of 0 to 0.8. The focus is on analyzing the effects of these variations on the total intensity ($G_{s}$) and the fundamental equilibrium solution.

Figure 2 visually illustrates the relationship between the total intensity $G_{s}$ and the depth of a uniform suspension, while incrementally increasing the forward scattering coefficient $A$ from 0 to 0.8 while keeping $\omega$ constant. As the value of $A$ increases, there is a noticeable increase in the intensity of $G_{s}$ at lower regions of the uniform suspension. Conversely, at the upper regions, the intensity of $G_{s}$ decreases.

In Figure 3(a), the plot illustrates the phototaxis function as a function of light intensity, considering a critical value of $G_{c}=1$. On the other hand, Figure 3(b) displays the influence of the forward scattering coefficient $A$ on the equilibrium position of the sublayer, using the same governing parameters.

In the case where $\omega=0$, the sublayer forms near the top of the suspension domain in the equilibrium state. However, as $A$ increases from 0 to 0.8, the sublayer’s equilibrium position gradually shifts towards the lower region.

Similarly, for $\omega=0.42$ and 0.47, the sublayer’s equilibrium position in the equilibrium state is located approximately at three-quarters and the middle height of the suspension, respectively. Once again, an increase in $A$ from 0 to 0.8 leads to a downward shift of the sublayer’s equilibrium position in each scenario.

To analyse stability, we employ linear perturbation theory, where a small perturbation with an amplitude $\epsilon\ll 1$ is introduced to the equilibrium state. This perturbation can be mathematically represented by the following equation:

$$[\bm{U},n,\zeta,L,<p>]=[0,n_{s},\zeta_{s},L_{s},<p_{s}>]+\epsilon[\bm{U}_{1},n_{1},\zeta_{1},L_{1},<\bm{p}_{1}>]+\mathcal{O}(\epsilon^{2})$$

(37)

By substituting the perturbed variables into equations (14) to (16) and linearizing the equations, we gather terms of $o(\epsilon)$ around the equilibrium state, leading to the following expression:

$$\bm{\nabla}\cdot\bm{U}_{1}=0,$$

(38)

where $\bm{U}_{1}=(U_{1},V_{1},W_{1})$.

$$Sc^{-1}\left(\frac{\partial\bm{U_{1}}}{\partial t}\right)+\sqrt{T_{a}}(z\times U_{1})+\bm{\nabla}P_{e}+Rn_{1}\hat{\bm{z}}=\nabla^{2}\bm{U_{1}},$$

(39)

$$\frac{\partial{n_{1}}}{\partial{t}}+V_{c}\bm{\nabla}\cdot(<\bm{p_{s}}>n_{1}+<\bm{p_{1}}>n_{s})+W_{1}\frac{dn_{s}}{dz}=\bm{\nabla}^{2}n_{1}.$$

(40)

If $G=G_{s}+\epsilon G_{1}+\mathcal{O}(\epsilon^{2})=(G_{s}^{c}+\epsilon G_{1}^{c})+(G_{s}^{d}+\epsilon G_{1}^{d})+\mathcal{O}(\epsilon^{2})$, then the steady collimated total intensity is perturbed as $L_{t}\exp\left(-\kappa\int_{z}^{1}(n_{s}(z^{\prime})+\epsilon n_{1}+\mathcal{O}(\epsilon^{2}))dz^{\prime}\right)$ and after simplification, we get

$$G_{1}^{c}=L_{t}\exp\left(-\int_{z}^{1}\kappa n_{s}(z^{\prime})dz^{\prime}\right)\left(\int_{1}^{z}\kappa n_{1}dz^{\prime}\right)$$

(41)

Similarly, $G_{1}^{d}$ is given by

$$G_{1}^{d}=\int_{0}^{4\pi}L_{1}^{d}(\bm{x},\bm{r})d\Omega.$$

(42)

Similarly, for the radiative heat flux $q=q_{s}+\epsilon q_{1}++\mathcal{O}(\epsilon^{2})=(q_{s}^{c}+\epsilon q_{1}^{c})+(q_{s}^{d}+\epsilon q_{1}^{d})+\mathcal{O}(\epsilon^{2})$, we find

$$\bm{q}_{1}^{c}=L_{t}\exp\left(-\int_{z}^{1}\kappa n_{s}(z^{\prime})dz^{\prime}\right)\left(\int_{1}^{z}\kappa n_{1}dz^{\prime}\right)\hat{z}$$

(43)

and

$$q_{1}^{d}=\int_{0}^{4\pi}L_{1}^{d}(\bm{x},\bm{r})\bm{r}d\Omega.$$

(44)

After perturbing the expression for swimming orientation and collecting $O(\epsilon)$ terms gives the perturbed swimming orientation

$$<\bm{p_{1}}>=G_{1}\frac{dM_{s}}{dG}\hat{\bm{z}}-M_{s}\frac{\bm{q_{1}}^{H}}{\bm{q_{s}}},$$

(45)

here, $\bm{q}_{1}^{H}=[\bm{q}_{1}^{x},\bm{q}_{1}^{y}]$ represents the horizontal component of the perturbed radiative flux $\bm{q}_{1}$. By substituting the value of $<\bm{p_{1}}>$ from Eq.$(\ref{46})$ into Eq.$(\ref{41})$ and simplifying, we obtain:

$$\frac{\partial{n_{1}}}{\partial{t}}+V_{c}\frac{\partial}{\partial z}\left(M_{s}n_{1}+n_{s}G_{1}\frac{dM_{s}}{dG}\right)-V_{c}n_{s}\frac{M_{s}}{q_{s}}\left(\frac{\partial q_{1}^{x}}{\partial x}+\frac{\partial q_{1}^{y}}{\partial y}\right)+W_{1}\frac{dn_{s}}{dz}=\nabla^{2}n_{1}.$$

(46)

By eliminating $P_{e}$ and the horizontal component of $u_{1}$, equations (39) and (40) can be simplified to three equations for the perturbed variables: the vertical component of the velocity $w_{1}$, the vertical component of the vorticity $\zeta_{1}$ (which is defined as $\zeta\cdot\hat{\bm{z}}$), and the concentration $n_{1}$. These variables can be further decomposed into normal modes as:

$$[W_{1},\zeta_{1},n_{1}]=[\tilde{W}(z),\tilde{Z}(z),\tilde{N}(z)]\exp{(\sigma t+i(lx+my))}.$$

(47)

The equation governing the perturbed intensity $L_{1}$ can be written as:

$$\xi\frac{\partial L_{1}}{\partial x}+\eta\frac{\partial L_{1}}{\partial y}+\nu\frac{\partial L_{1}}{\partial z}+\kappa(n_{s}L_{1}+n_{1}L_{s})=\frac{\omega\kappa}{4\pi}(n_{s}M_{1}+G_{s}n_{1}+A\nu(n_{s}q_{s}\cdot\hat{z}-q_{s}n_{1}))-\kappa L_{s}n_{1},$$

(48)

with boundary conditions

	$$L_{1}(x,y,z=1,\xi,\eta,\nu)=0,\qquad\theta\in[\pi/2,\pi],~{}~{}\phi\in[0,2\pi],$$		(49a)
	$$L_{1}(x,y,z=0,\xi,\eta,\nu)=0,\qquad\theta\in[0,\pi/2],~{}~{}\phi\in[0,2\pi].$$		(49b)

The Eq. $(\ref{48})$ implies that $L_{1}^{d}$ can be represented by the following expression:

$$L_{1}^{d}=\Psi^{d}(z,\xi,\eta,\nu)\exp{(\sigma t+i(lx+my))}.$$

From Eqs. (41) and (42), we get

$$G_{1}^{c}=\left[L_{t}\exp\left(-\int_{z}^{1}\kappa n_{s}(z^{\prime})dz^{\prime}\right)\left(\int_{1}^{z}\kappa n_{1}dz^{\prime}\right)\right]\exp{(\sigma t+i(lx+my))}=\mathbb{G}^{c}(z)\exp{(\sigma t+i(lx+my))},$$

(50)

and

$$G_{1}^{d}=\mathbb{G}^{d}(z)\exp{(\sigma t+i(lx+my))}=\left(\int_{0}^{4\pi}\Psi^{d}(z,\xi,\eta,\nu)d\Omega\right)\exp{(\sigma t+i(lx+my))},$$

(51)

where $\mathbb{G}(z)=\mathbb{G}^{c}(z)+\mathbb{G}^{d}(z)$ is the perturbed total intensity. Similarly from Eqs. (43) and (44), we have

$$\bm{q}_{1}=[q_{1}^{x},q_{1}^{y},q_{1}^{z}]=[P(z),Q(z),S(z)]\exp{[\sigma t+i(lx+my)]},$$

where

$$[P(z),Q(z)]=\int_{0}^{4\pi}[\xi,\eta]\Psi^{d}(z,\xi,\eta,\nu)d\Omega,$$

and

$$S(z)=\int_{0}^{4\pi}[\Psi^{c}(z,\xi,\eta,\nu)+\Psi^{d}(z,\xi,\eta,\nu)]\nu d\Omega.$$

Note that the $P(z)$ and $Q(z)$ appear due to scattering. On the other hand, $S(z)$ appears due to anisotropic scattering which becomes zero in the case of isotropic scattering.

Now $\Psi^{d}$ satisfies

$$\frac{d\Psi^{d}}{dz}+\frac{(i(l\xi+m\eta)+\kappa n_{s})}{\nu}\Psi^{d}=\frac{\omega\kappa}{4\pi\nu}(n_{s}\mathcal{\mathbb{G}}+G_{s}\tilde{N}(z)+A\nu(n_{s}S-q_{s}\tilde{N}(z)))-\frac{\kappa}{\nu}I_{s}\tilde{N}(z),$$

(52)

subject to the boundary conditions

	$$\Psi^{d}(1,\xi,\eta,\nu)=0,\qquad\theta\in[\pi/2,\pi],~{}~{}\phi\in[0,2\pi],$$		(53a)
	$$\Psi^{d}(0,\xi,\eta,\nu)=0,\qquad\theta\in[0,\pi/2],~{}~{}\phi\in[0,2\pi].$$		(53b)

The stability equations reformed as

$$\left(\sigma S_{c}^{-1}+k^{2}-\frac{d^{2}}{dz^{2}}\right)\left(\frac{d^{2}}{dz^{2}}-k^{2}\right)\tilde{W}=Rk^{2}\tilde{N}(z),$$

(54)

$$\left(\sigma S_{c}^{-1}+k^{2}-\frac{d^{2}}{dz^{2}}\right)\tilde{Z}(z)=\sqrt{T_{a}}\frac{d\tilde{W}}{dz}$$

(55)

$$\left(\sigma+k^{2}-\frac{d^{2}}{dz^{2}}\right)\tilde{N}(z)+V_{c}\frac{d}{dz}\left(M_{s}\tilde{N}(z)+n_{s}\mathbb{G}\frac{dM_{s}}{dG}\right)-i\frac{V_{c}n_{s}M_{s}}{q_{s}}(lP+mQ)=-\frac{dn_{s}}{dz}\tilde{W}(z),$$

(56)

with boundary conditions

$$\tilde{W}(z)=\frac{d\tilde{W}(z)}{dz}=\tilde{Z}(z)=\frac{d\tilde{N}(z)}{dz}-V_{c}M_{s}\tilde{N}(z)-n_{s}V_{c}\mathbb{G}\frac{dM_{s}}{dG}=0\quad at\quad z=0,$$

(57)

$$\tilde{W}(z)=\frac{d\tilde{W}(z)}{dz}=\frac{d\tilde{Z}(z)}{dz}=\frac{d\tilde{N}(z)}{dz}-V_{c}M_{s}\tilde{N}(z)-n_{s}V_{c}\mathbb{G}\frac{dM_{s}}{dG}=0\quad at\quad z=1.$$

(58)

n the given context, the non-dimensional wavenumber is represented by $k$, which is determined by taking the square root of the sum of the squares of $l$ and $m$. Equations (57)-(58) form an eigenvalue problem, where $\sigma$ is described as a function of several dimensionless parameters such as $V_{c}$, $\kappa$, $\omega$, $A$, $l$, $m$, and $R$. Equation (56) can be equivalently expressed as follows:

$$D^{2}\tilde{N}(z)-\Lambda_{3}(z)D\tilde{N}(z)-(\sigma+k^{2}+\Lambda_{2}(z))\tilde{N}(z)-\Lambda_{1}(z)\int_{1}^{z}\tilde{N}(z)dz-\Lambda_{0}(z)=Dn_{s}\tilde{W},$$

(59)

where

	$$\Lambda_{0}(z)=V_{c}D\left(n_{s}\mathbb{G}^{d}\frac{dM_{s}}{dG}\right)-\iota\frac{V_{c}n_{s}M_{s}}{q_{s}}(lP+mQ),$$		(60a)
	$$\Lambda(z)=\kappa V_{c}D\left(n_{s}G_{s}^{c}\frac{dM_{s}}{dG}\right)$$		(60b)
	$$\Lambda_{2}(z)=2\kappa V_{c}n_{s}G_{s}^{c}\frac{dM_{s}}{dG}+V_{c}\frac{dM_{s}}{dM}DG_{s}^{d},$$		(60c)
	$$\Lambda_{3}(z)=V_{c}M_{s}.$$		(60d)

Introducing a novel variable denoted as

$$\tilde{\Theta}(z)=\int_{1}^{z}\tilde{N}(z^{\prime})dz^{\prime},$$

(61)

the linear stability equations can be reformulated as

$$D^{4}\tilde{W}-\left(\sigma S_{c}^{-1}+k^{2}\right)D^{2}\tilde{W}-\left(\sigma S_{c}^{-1}+k^{2}\right)\tilde{W}=Rk^{2}D\tilde{\Theta},$$

(62)

$$D^{2}\tilde{Z}(z)-\left(\sigma S_{c}^{-1}+k^{2}-D^{2}\right)\tilde{Z}(z)=\sqrt{T_{a}}D\tilde{W}$$

(63)

$$D^{3}\tilde{\Theta}-\Lambda_{3}(z)D^{2}\tilde{\Theta}-(\sigma+k^{2}+\Lambda_{2}(z))D\tilde{\Theta}-\Lambda_{1}(z)\tilde{\Theta}-\Lambda_{0}(z)=Dn_{s}\tilde{W}.$$

(64)

The boundary conditions remain unchanged except for the flux boundary condition, given by the equation

$$D^{2}\tilde{N}(z)-\Lambda_{2}D\tilde{N}(z)-n_{s}V_{c}\mathbb{G}\frac{dM_{s}}{dG}=0\quad at\quad z=0,1.$$

(65)

Furthermore, an additional boundary condition is introduced

$$\tilde{\Theta}(z)=0,\quad at\quad z=1.$$

(66)

To obtain the neutral (marginal) stability curves or the growth rate $Re(\sigma)$ as a function of R in the (k, R)-plane for a fixed parameter set, a fourth-order accurate finite-difference scheme based on Newton-Raphson-Kantorovich (NRK) iterations is utilized to solve equations 62 to 64. By analysing the resulting graph, points where $Re(\sigma)=0$ are identified, representing the marginal (neutral) stability curve. On this curve, if $Im(\sigma)=0$, it indicates a stationary (non-oscillatory) bioconvective solution, while non-zero $Im(\sigma)$ indicates oscillatory solutions. If the most unstable mode remains on the oscillatory branch of the neutral curve, it signifies overstability.

When oscillatory solutions are present, there is a common point $(k_{b})$ where the stationary and oscillatory branches intersect, and the oscillatory branch includes points with $k\leq k_{b}$. Among the points on the neutral curve $R^{(n)}(k)$ (where n = 1, 2, 3, …), a particular most unstable mode $(k_{c},R_{c})$ is selected. The wavelength of the initial disturbance is then determined as $\lambda_{c}=2\pi/k_{c}$. The bioconvective solution is referred to as mode n if it can be organized into n convection cells such that one cell overlaps another vertically.

To identify the most unstable mode starting from an initial equilibrium solution, we employ a specific set of fixed parameters. These parameters include $S_{c}=20$, $G_{c}=1$, $V_{c}=10$, $15$, $20$, with $\omega$ ranging from $0$ to $1$, $\kappa$ values of $0.5$ and $1$, and $A$ values of $0$, $0.4$, and $0.8$. The Taylor number, which represents the rotation rate, is varied from lower non-zero values to higher values. To investigate the impact of rotation on a forward scattering algal suspension, we focus on specific values of the scattering albedo $\omega$. For the case of $\kappa=0.5$, we examine the scenarios where $\omega$ takes on the values $0.1$, $0.42$, and $0.47$. Similarly, for the case of $\kappa=1$, we consider $\omega$ values of $0.1$, $0.58$, and $0.60$. When $\omega=0.1$, the equilibrium state’s sublayer is positioned near the top of the domain. As the value of $\omega$ increases, the sublayer’s equilibrium position shifts from the top to approximately three-quarters height, and then further to the middle height of the domain.