2.1 Governing equations

Owing to the small scales relevant to the swimming of micro-organisms, the Reynolds number is small and the inertia of the fluid is not important. We thus neglect inertia in our analysis. We now begin by writing the mass conservation and the momentum conservation equations in the fluid as follows:

(2.2)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}y}=0, & \displaystyle\end{eqnarray}$$

(2.3)$$\begin{eqnarray}\displaystyle & \displaystyle 0=-\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left(\unicode[STIX]{x1D707}\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}x}\right)+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}\left(\unicode[STIX]{x1D707}\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}y}\right), & \displaystyle\end{eqnarray}$$

(2.4)$$\begin{eqnarray}\displaystyle & \displaystyle 0=-\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}y}-\unicode[STIX]{x1D70C}_{\infty }g+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left(\unicode[STIX]{x1D707}\frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}x}\right)+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}\left(\unicode[STIX]{x1D707}\frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}y}\right). & \displaystyle\end{eqnarray}$$

Here, $u$ and $v$ represent the velocity fields in the $x$ and $y$ directions, respectively. The equations are written in the frame of reference of the swimming sheet, which is translating with a constant velocity $V$ in the positive $y$ direction.

In addition, we have the advection–diffusion equation describing the distribution of the stratifying agent. For a linear variation of viscosity with the temperature/concentration, the advection–diffusion equation can be written in terms of the fluid viscosity as follows:

(2.5)$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}(\unicode[STIX]{x1D707}u)+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}(\unicode[STIX]{x1D707}v)=c\left(\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D707}}{\unicode[STIX]{x2202}x^{2}}+\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D707}}{\unicode[STIX]{x2202}y^{2}}\right).\end{eqnarray}$$

Similar equations are found in the literature for studies involving the transport of particles in a density-stratified fluid, where the density varies linearly with the temperature/concentration (Candelier, Mehaddi & Vauquelin Reference Candelier, Mehaddi and Vauquelin2014). Here, $c$ is the diffusion coefficient. The viscosity and the pressure fields can be written as a combination of undisturbed and disturbed fields as $\unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}_{0}+\unicode[STIX]{x1D707}^{\prime }$ and $p=p_{0}+p^{\prime }$, respectively. We note that, according to the Boussinesq approximation, for small changes in the temperature/concentration fields, the material properties can be approximated by their values evaluated at an ambient temperature/concentration except for the body-force term (Leal Reference Leal2007). Hence, under the Boussinesq approximation, the viscosity field in the momentum conservation equations can be approximated by the reference fluid viscosity $\unicode[STIX]{x1D707}_{\infty }$. However, here in this work, in addition to the reference fluid viscosity, we retain the disturbed viscosity field due to the swimmer’s motion ($\unicode[STIX]{x1D707}^{\prime }$) in our analysis, as we expect the swimmer dynamics to be influenced by the disturbed viscosity field. Consequently, we represent the viscosity field as $\unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}_{\infty }+\unicode[STIX]{x1D707}^{\prime }$.

Next, we substitute these in the governing equations of motion and make the equations dimensionless using the following scales: length scale $l_{c}=1/k$, velocity scale $u_{c}=\unicode[STIX]{x1D70E}/k$, pressure scale $p_{c}=\unicode[STIX]{x1D707}_{\infty }\unicode[STIX]{x1D70E}$ and time scale $t_{c}=1/\unicode[STIX]{x1D70E}$. Note that the undisturbed pressure gradient $\unicode[STIX]{x2202}p_{0}/\unicode[STIX]{x2202}y$ will be balanced by the body-force term $\unicode[STIX]{x1D70C}_{\infty }g$. Also, in the frame of reference of the swimmer, the undisturbed fluid viscosity is changing spatially as well as temporally, i.e. $\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{0}/\unicode[STIX]{x2202}y=-\unicode[STIX]{x1D6FE}$ and $\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{0}/\unicode[STIX]{x2202}t=\unicode[STIX]{x1D6FE}V$.

We also use the streamfunction variable to represent the velocity fields as follows: $u=-\unicode[STIX]{x2202}\unicode[STIX]{x1D713}/\unicode[STIX]{x2202}y$ and $v=\unicode[STIX]{x2202}\unicode[STIX]{x1D713}/\unicode[STIX]{x2202}x$. In addition, we eliminate the pressure terms in (2.3) and (2.4) to obtain a single equation in terms of the streamfunction. We use a variable transformation $z=y+t$ to simplify our calculations. Note that such a transformation has been previously performed in the literature for problems involving the motion of a swimming sheet (Lauga Reference Lauga2007). After performing the above operations, we obtain the following equations of motion:

(2.6)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FB}^{4}\unicode[STIX]{x1D713} & = & \displaystyle -\unicode[STIX]{x1D6FC}\left(\unicode[STIX]{x1D707}^{\prime }\unicode[STIX]{x1D6FB}^{4}\unicode[STIX]{x1D713}+2\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}^{\prime }}{\unicode[STIX]{x2202}x}\left(\frac{\unicode[STIX]{x2202}^{3}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}x^{3}}+\frac{\unicode[STIX]{x2202}^{3}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}x\unicode[STIX]{x2202}z^{2}}\right)+2\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}^{\prime }}{\unicode[STIX]{x2202}z}\left(\frac{\unicode[STIX]{x2202}^{3}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}z^{3}}+\frac{\unicode[STIX]{x2202}^{3}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}z\unicode[STIX]{x2202}x^{2}}\right)\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D707}^{\prime }}{\unicode[STIX]{x2202}x^{2}}\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}x^{2}}+\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D707}^{\prime }}{\unicode[STIX]{x2202}z^{2}}\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}z^{2}}+2\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D707}^{\prime }}{\unicode[STIX]{x2202}x\unicode[STIX]{x2202}z}\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}x\unicode[STIX]{x2202}z}\right).\end{eqnarray}$$

(2.7)$$\begin{eqnarray}Pe\left(\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}^{\prime }}{\unicode[STIX]{x2202}z}-V\right)+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left(-\unicode[STIX]{x1D707}^{\prime }\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}z}\right)+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}\left(\unicode[STIX]{x1D707}^{\prime }\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}x}\right)-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}x}\right)=\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D707}^{\prime }}{\unicode[STIX]{x2202}x^{2}}+\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D707}^{\prime }}{\unicode[STIX]{x2202}z^{2}}.\end{eqnarray}$$

Here, $Pe=\unicode[STIX]{x1D70E}/k^{2}c$ is the Péclet number, which signifies the ratio between convection and diffusion of the stratifying agent, and $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D6FE}/\unicode[STIX]{x1D707}_{\infty }k$ is a dimensionless parameter, which signifies the strength of the viscosity stratification. We aim to solve the above system of two coupled partial differential equations with two unknown variables ($\unicode[STIX]{x1D713},\unicode[STIX]{x1D707}^{\prime }$).

2.2 Boundary conditions

On the sheet surface, we apply a no-slip boundary condition, which is given by

(2.8a,b)$$\begin{eqnarray}u|_{x=bk\sin (z)}=bk\cos (z),\quad v|_{x=bk\sin (z)}=0.\end{eqnarray}$$

In a stationary frame of reference, the fluid velocity is zero far away from the sheet. Hence, in the frame of reference of the sheet, we have

(2.9a,b)$$\begin{eqnarray}u|_{x\rightarrow \infty }=0,\quad v|_{x\rightarrow \infty }=-V.\end{eqnarray}$$

Additionally, we apply a no-flux boundary condition at the surface of the swimmer, which implies that the swimmer surface is impermeable when concentration gradients are present:

(2.10)$$\begin{eqnarray}\unicode[STIX]{x1D735}(\unicode[STIX]{x1D707}_{0}+\unicode[STIX]{x1D707}^{\prime })\boldsymbol{\cdot }\boldsymbol{n}|_{x=bk\sin (z)}=0.\end{eqnarray}$$

Here, $\boldsymbol{n}$ is the normal vector calculated at the surface of the swimmer.

As the disturbed viscosity should vanish far away from the swimmer, we have the following:

(2.11)$$\begin{eqnarray}\unicode[STIX]{x1D707}^{\prime }|_{x\rightarrow \infty }=0.\end{eqnarray}$$

3.1 Solution at $O(bk)$

The first-order governing equations are given as follows:

(3.1a)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}^{4}\unicode[STIX]{x1D713}_{1}=0, & \displaystyle\end{eqnarray}$$

(3.1b)$$\begin{eqnarray}\displaystyle & \displaystyle Pe\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{1}^{\prime }}{\unicode[STIX]{x2202}x}-V_{1}-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}x}\right)=\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D707}_{1}^{\prime }}{\unicode[STIX]{x2202}x^{2}}+\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D707}_{1}^{\prime }}{\unicode[STIX]{x2202}z^{2}}. & \displaystyle\end{eqnarray}$$

$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}x}=0,\quad \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}z}=-\cos (z),\quad \text{at }x=0, & \displaystyle \nonumber\\ \displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{1}^{\prime }}{\unicode[STIX]{x2202}x}=-\cos (z),\quad \text{at }x=0, & \displaystyle \nonumber\\ \displaystyle & \displaystyle \unicode[STIX]{x1D707}_{1}^{\prime }\rightarrow 0,\quad \text{as }x\rightarrow \infty . & \displaystyle \nonumber\end{eqnarray}$$

We observe that the boundary conditions are periodic functions in the vertical direction with a mean value of zero. Hence, we expect the first-order streamfunction and viscosity to be periodic functions as well with a zero mean in the vertical direction. Now, as all the terms in (3.1b) have a zero vertical mean, we expect the first-order velocity to be zero as well. This implies the first-order sheet velocity is zero in a viscosity-stratified fluid, which is the same as that in a homogeneous fluid (Taylor Reference Taylor1951).

We now determine the solutions for $\unicode[STIX]{x1D713}_{1}$ and $\unicode[STIX]{x1D707}_{1}^{\prime }$ by solving (3.1) and satisfying the boundary conditions. The first-order streamfunction, which is governed by the biharmonic equation, can be found to be as follows:

(3.2)$$\begin{eqnarray}\unicode[STIX]{x1D713}_{1}=-(1+x)\sin (z)\text{e}^{-x}.\end{eqnarray}$$

We substitute this in (3.1b) to obtain the first-order disturbed viscosity as follows:

(3.3)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D707}_{1}^{\prime } & = & \displaystyle \text{Re}\left(\frac{\text{i}}{Pe\sqrt{1+\text{i}Pe}}\text{e}^{-\sqrt{1+\text{i}Pe}\,x}(\cos (z)+\text{i}\sin (z))\right.\nonumber\\ \displaystyle & & \displaystyle -\left.\frac{\text{i}}{Pe\sqrt{1-\text{i}Pe}}\text{e}^{-\sqrt{1-\text{i}Pe}\,x}(\cos (z)-\text{i}\sin (z))-x\text{e}^{-x}\cos (z)+\frac{2}{Pe}\text{e}^{-x}\sin (z)\right).\qquad\end{eqnarray}$$

Here, $\text{i}^{2}=-1$ and $\text{Re}(\cdot )$ denotes the real part of a function.

3.2 Solution at $O(bk)^{2}$

The second-order governing equations are given as follows:

(3.4a)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FB}^{4}\unicode[STIX]{x1D713}_{2} & = & \displaystyle -\unicode[STIX]{x1D6FC}\left(\unicode[STIX]{x1D707}_{1}^{\prime }\unicode[STIX]{x1D6FB}^{4}\unicode[STIX]{x1D713}_{1}+2\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{1}^{\prime }}{\unicode[STIX]{x2202}x}\left(\frac{\unicode[STIX]{x2202}^{3}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}x^{3}}+\frac{\unicode[STIX]{x2202}^{3}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}x\unicode[STIX]{x2202}z^{2}}\right)+2\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{1}^{\prime }}{\unicode[STIX]{x2202}z}\left(\frac{\unicode[STIX]{x2202}^{3}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}z^{3}}+\frac{\unicode[STIX]{x2202}^{3}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}z\unicode[STIX]{x2202}x^{2}}\right)\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D707}_{1}^{\prime }}{\unicode[STIX]{x2202}x^{2}}\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}x^{2}}+\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D707}_{1}^{\prime }}{\unicode[STIX]{x2202}z^{2}}\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}z^{2}}+2\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D707}_{1}^{\prime }}{\unicode[STIX]{x2202}x\unicode[STIX]{x2202}z}\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}x\unicode[STIX]{x2202}z}\right),\end{eqnarray}$$

(3.4b)$$\begin{eqnarray}Pe\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{2}^{\prime }}{\unicode[STIX]{x2202}z}-V_{2}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left(-\unicode[STIX]{x1D707}_{1}^{\prime }\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}z}\right)+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}\left(\unicode[STIX]{x1D707}_{1}^{\prime }\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}x}\right)-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{2}}{\unicode[STIX]{x2202}x}\right)=\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D707}_{2}^{\prime }}{\unicode[STIX]{x2202}x^{2}}+\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D707}_{2}^{\prime }}{\unicode[STIX]{x2202}z^{2}}.\end{eqnarray}$$

(3.5a,b)$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{2}}{\unicode[STIX]{x2202}x}=-\text{sin}(z)\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}x^{2}},\quad \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{2}}{\unicode[STIX]{x2202}z}=-\text{sin}(z)\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}x\unicode[STIX]{x2202}z},\quad \text{at }x=0.\end{eqnarray}$$

(3.6)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{2}^{\prime }}{\unicode[STIX]{x2202}x}=-\text{sin}(z)\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D707}_{1}^{\prime }}{\unicode[STIX]{x2202}x^{2}}+\cos (z)\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{1}^{\prime }}{\unicode[STIX]{x2202}z},\quad \text{at }x=0. & \displaystyle\end{eqnarray}$$

(3.7)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D707}_{2}^{\prime }\rightarrow 0,\quad \text{as }x\rightarrow \infty . & \displaystyle\end{eqnarray}$$

We find that the vertical mean value of second-order vertical velocity and the viscosity gradient evaluated at the surface of the sheet are non-zero. Consequently, the velocity and viscosity fields must have a non-zero mean component to satisfy the boundary conditions at the sheet surface. To obtain these mean components, we define the vertical averaged viscosity and streamfunction in the domain as follows:

(3.8a,b)$$\begin{eqnarray}\langle \unicode[STIX]{x1D713}_{2}\rangle =\frac{1}{2\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}\unicode[STIX]{x1D713}_{2}(x,z)\,\text{d}z,\quad \langle \unicode[STIX]{x1D707}_{2}^{\prime }\rangle =\frac{1}{2\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}\unicode[STIX]{x1D707}_{2}^{\prime }(x,z)\,\text{d}z.\end{eqnarray}$$

We now average over the vertical direction on both sides of (3.4) to obtain

(3.9a)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}^{4}\langle \unicode[STIX]{x1D713}_{2}\rangle }{\text{d}x^{4}}=-\unicode[STIX]{x1D6FC}k_{1}(x), & \displaystyle\end{eqnarray}$$

(3.9b)$$\begin{eqnarray}\displaystyle & \displaystyle Pe\left(-V_{2}+k_{2}(x)-\frac{\text{d}\langle \unicode[STIX]{x1D713}_{2}\rangle }{\text{d}x}\right)=\frac{\text{d}^{2}\langle \unicode[STIX]{x1D707}_{2}^{\prime }\rangle }{\text{d}x^{2}}, & \displaystyle\end{eqnarray}$$

$$\begin{eqnarray}\displaystyle k_{1}(x) & = & \displaystyle \left\langle \unicode[STIX]{x1D707}_{1}^{\prime }\unicode[STIX]{x1D6FB}^{4}\unicode[STIX]{x1D713}_{1}+2\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{1}^{\prime }}{\unicode[STIX]{x2202}x}\left(\frac{\unicode[STIX]{x2202}^{3}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}x^{3}}+\frac{\unicode[STIX]{x2202}^{3}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}x\unicode[STIX]{x2202}z^{2}}\right)+2\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{1}^{\prime }}{\unicode[STIX]{x2202}z}\left(\frac{\unicode[STIX]{x2202}^{3}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}z^{3}}+\frac{\unicode[STIX]{x2202}^{3}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}z\unicode[STIX]{x2202}x^{2}}\right)+\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D707}_{1}^{\prime }}{\unicode[STIX]{x2202}x^{2}}\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}x^{2}}\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D707}_{1}^{\prime }}{\unicode[STIX]{x2202}z^{2}}\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}z^{2}}+2\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D707}_{1}^{\prime }}{\unicode[STIX]{x2202}x\unicode[STIX]{x2202}z}\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}x\unicode[STIX]{x2202}z}\right\rangle \nonumber\end{eqnarray}$$

and

$$\begin{eqnarray}k_{2}(x)=\left\langle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left(-\unicode[STIX]{x1D707}_{1}^{\prime }\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}z}\right)+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}\left(\unicode[STIX]{x1D707}_{1}^{\prime }\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}x}\right)\right\rangle .\end{eqnarray}$$

These expressions can be obtained as we have determined the first-order velocity and viscosity fields.

We now define $\langle \unicode[STIX]{x1D713}_{2}\rangle =-V_{2}x+\unicode[STIX]{x1D713}_{2s}$, where $\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{2s}/\unicode[STIX]{x2202}x$ is the $y$ direction velocity field, averaged in the vertical direction and measured in a stationary frame of reference. We now use the vertical velocity boundary condition (refer to (3.5)) to obtain the expression for the swimming sheet velocity as follows:

(3.10)$$\begin{eqnarray}V_{2}=\left\langle \sin (z)\frac{\text{d}^{2}\unicode[STIX]{x1D713}_{1}}{\text{d}x^{2}}(x=0)\right\rangle +\frac{\text{d}\unicode[STIX]{x1D713}_{2s}}{\text{d}x}(x=0).\end{eqnarray}$$

Here, $\unicode[STIX]{x1D713}_{1}$ and $\unicode[STIX]{x1D713}_{2s}$ are determined by solving (3.2) and (3.9a), respectively. After substituting these expressions, we obtain the velocity of a swimmer in a viscosity-stratified fluid as

(3.11)$$\begin{eqnarray}V_{2}=\frac{\unicode[STIX]{x1D70E}}{k}b^{2}k^{2}\left(\frac{1}{2}+\unicode[STIX]{x1D6FC}\,\text{Re}[f(Pe)]\right),\end{eqnarray}$$

where the function $f(Pe)$ is given by

(3.12)$$\begin{eqnarray}\displaystyle f(Pe) & = & \displaystyle \frac{1}{8(1+\sqrt{1-\text{i}Pe})^{2}(1+\sqrt{1+\text{i}Pe})^{2}(1+Pe^{2})Pe^{3}}\nonumber\\ \displaystyle & & \displaystyle \times \{(112\text{i}Pe+14\text{i}Pe^{3}-2\text{i}Pe^{5}-44Pe^{2})(\sqrt{1-\text{i}Pe}-\sqrt{1+\text{i}Pe})\nonumber\\ \displaystyle & & \displaystyle -\,Pe^{6}-4Pe^{2}(49-31\sqrt{1+Pe^{2}})\nonumber\\ \displaystyle & & \displaystyle +\,Pe^{4}(-5-4\sqrt{1+Pe^{2}})+192(-1+\sqrt{1+Pe^{2}})\nonumber\\ \displaystyle & & \displaystyle +\,(\sqrt{1-\text{i}Pe}+\sqrt{1+\text{i}Pe})[-4Pe^{2}+4Pe^{4}-48(-1+\sqrt{1+Pe^{2}})]\!\}.\qquad\end{eqnarray}$$

We expand the above equation using a Taylor series expansion approach for the cases of both low and high Péclet numbers to obtain the following.

Low Péclet numbers $(Pe\ll 1)$:

(3.13)$$\begin{eqnarray}V_{2}=\frac{\unicode[STIX]{x1D70E}}{k}b^{2}k^{2}\left(\frac{1}{2}+\unicode[STIX]{x1D6FC}\left(\frac{Pe}{4}-\frac{741}{2048}Pe^{3}+\frac{3289}{8192}Pe^{5}+O(Pe^{6})\right)\right).\end{eqnarray}$$

High Péclet numbers $(Pe\gg 1)$:

(3.14)$$\begin{eqnarray}V_{2}=\frac{\unicode[STIX]{x1D70E}}{k}b^{2}k^{2}\left(\frac{1}{2}+\unicode[STIX]{x1D6FC}\left(-\frac{1}{8Pe}+O(Pe^{-5/2})\right)\right).\end{eqnarray}$$

4.1 Effect of viscosity gradient on the swimming velocity

The velocity of the swimming sheet in a fluid with a linear viscosity profile is given in (3.11). In addition to the homogeneous contribution to the swimming velocity, the contribution of the viscosity gradient to the non-dimensional velocity of the swimmer can be expressed as

(4.1)$$\begin{eqnarray}V_{viscous}=\unicode[STIX]{x1D6FC}\,\text{Re}[f(Pe)].\end{eqnarray}$$

From the above equation, we find that the velocity is governed by two parameters, $\unicode[STIX]{x1D6FC}$, which signifies the strength of the viscosity gradient, and $Pe$, which signifies the relative strength of the advection and diffusion of the stratifying agent. The dependence on $\unicode[STIX]{x1D6FC}$ is linear, which implies that a stronger viscosity stratification would enhance the contribution of $V_{viscous}$ to the velocity of the swimmer. The dependence on $Pe$ is given by (3.12), which is plotted in figure 2. From figure 2, we find that the linear fit derived from (3.13) ($f(Pe)=Pe/4$), provides an accurate estimate for Péclet numbers up to 0.35. We also observe that $|\,f(Pe)|$ has a maximum value of ${\sim}0.1$. Additionally, as the spatial scale on which the magnitude of the stratifying agent varies ($\unicode[STIX]{x1D707}_{\infty }/\unicode[STIX]{x1D6FE}$) is expected to be much larger than the length scale of the swimmer ($1/k$), we expect the value of $\unicode[STIX]{x1D6FC}$ to be much smaller than unity. Hence, the contribution of $V_{viscous}$ to the swimming velocity (refer to (4.1)) is expected to be at least two orders of magnitude smaller compared to the homogeneous counterpart (refer to (3.11)). Nevertheless, by focusing on the contribution of the viscosity-gradient-induced velocity, we can obtain important insights regarding the preferential bias of the swimmers towards high/low-viscosity regions.

We observe that, for small Péclet numbers ($0<Pe<3$), the effect of the viscosity gradient is to propel the swimmer in the positive $y$ direction. However, for large Péclet numbers ($Pe>3$), the viscosity gradient pushes the swimmer in the negative $y$ direction. The linear dependence of $V_{viscous}$ on $\unicode[STIX]{x1D6FC}$ implies that if the direction of the viscosity gradient is reversed, the effect of the viscosity gradient on the swimming velocity would be reversed. Hence, for small (large) Péclet numbers, the effect of the viscosity gradient is to propel the swimmer towards regions of low (high) viscosity.

Figure 2.

Contribution of the viscosity gradient to the swimming velocity for $\unicode[STIX]{x1D6FC}=1$. The blue line represents the linear fit at small Péclet numbers derived from (3.13).

In our analysis, we have assumed a surface transverse wave travelling in the negative $y$ direction. However, we can show that for a wave travelling in the positive $y$ direction, the expression for $V_{viscous}$ remains the same. This implies that the contribution of the viscosity gradient to the swimming velocity is independent of the direction of the travelling wave generated by the swimmer. We can now conclude that, in the presence of viscosity gradients, the swimming sheet would prefer regions of low (high) viscosity for small (large) Péclet numbers. Consequently, our analysis shows that the swimming sheet shows negative viscotaxis for $0<Pe<3$ and positive viscotaxis for $Pe>3$. As opposed to the previous studies, which attributed the viscotaxis response to the speculative existence of viscoreceptors (Petrino & Doetsch Reference Petrino and Doetsch1978; Daniels et al. Reference Daniels, Longland and Gilbart1980; Sherman et al. Reference Sherman, Timkina and Glagolev1982; Takabe et al. Reference Takabe, Tahara, Islam, Affroze, Kudo and Nakamura2017), our work shows that swimmers might exhibit a viscotaxis behaviour purely based on hydrodynamic effects.

The hydrodynamic mechanism for such a behaviour can be explained by focusing on the distribution of the viscosity in the vicinity of the swimmer. Figure 3(a,b) shows the variation of the disturbed viscosity near the sheet’s crest (shown in blue) for the cases of low Péclet number ($Pe=0.01$) and high Péclet number ($Pe=100$), respectively. For low Péclet numbers, we find that the value of the disturbed viscosity ($\unicode[STIX]{x1D707}_{1}$) is positive in the lower half of the sheet surface ($0<y<\unicode[STIX]{x03C0}/2$) and negative in the upper half ($\unicode[STIX]{x03C0}/2<y<\unicode[STIX]{x03C0}$). However, for large Péclet numbers, this trend is reversed and even the magnitude of the viscosity is smaller. As the sheet retreats back towards its undeformed configuration, viscous forces are generated in the fluid which resist the velocity gradients. The direction of these viscous forces is shown by the green arrows in figure 3(c,d). For low (high) Péclet numbers, the lower (upper) half of the sheet’s crest will experience more viscous resistance due to the higher fluid viscosity. These regions are shown by the longer green arrows in figure 3(c,d). This asymmetric distribution of the viscous forces creates a vortical fluid flow around the sheet’s crest which propels the sheet (red arrows indicate the direction of propulsion). We note that the direction of the vortical flow reverses at high Péclet numbers, which causes the contribution of the viscosity gradient to the swimming velocity to be negative for high $Pe$ (refer to figure 2).

Figure 3.

Distribution of the disturbed viscosity ($\unicode[STIX]{x1D707}_{1}$) for (a) $Pe=0.01$ and (b) $Pe=100$. Hydrodynamic mechanism describing the contribution of viscosity gradients to the sheet’s motion for (c) $Pe=0.01$ and (d) $Pe=100$. The ‘$+$’ sign indicates a positive value for the disturbed viscosity and the ‘$-$’ sign indicates a negative value for the disturbed viscosity; green arrows represent the viscous forces in the fluid; and red arrows represent the contribution of the vortical flow (shown in green) to the propulsion of the swimmer.

Studies have shown that diffusive transport is significantly more important compared to advective transport within the mucous layer (for example, Mirbagheri & Fu (Reference Mirbagheri and Fu2016) reported a Péclet number of $10^{-3}$ for transport of ammonia in the gastric mucous layer, and Lewis, Keener & Fogelson (Reference Lewis, Keener and Fogelson2017) reported a Péclet number of $10^{-2}$ for transport of hydrogen ions within the gastric mucous layer). Additionally, diffusion is the predominant transport process in biofilm cells, where the Péclet number is of the order of $10^{-3}$ (Stewart Reference Stewart2012). For such small Péclet numbers observed in mucous layers or biofilms, our analysis predicts that the swimmers would prefer regions of low viscosities.

Additionally, when our analysis is extended for large viscosity gradients, the contribution of the viscosity gradient dominates the homogeneous contribution and our theory predicts accumulation of the swimmers in regions of low or high viscosity depending on the Péclet number corresponding to the transport of the stratifying agent. Stehnach (Reference Stehnach2019) observed that swimming biflagellates accumulate in the high-viscosity region for longer flagella lengths. For their experimental analysis, the length scale is $l_{c}\sim 10^{-5}~\text{m}$, the velocity scale is $u_{c}\sim 10^{-4}~\text{m}~\text{s}^{-1}$ and the diffusion coefficient is $c\sim 10^{-11}~\text{m}^{2}~\text{s}^{-1}$, which results in a Péclet number of $Pe\sim 100$. If we assume that the swimming sheet model can be used to provide an insight into the motion of C. reinhardtii in a viscosity-stratified fluid, then from our analysis we can predict that, for $Pe\sim 100$, the swimmers would accumulate in the high-viscosity regions as observed in their experiments.

4.2 Power expenditure and swimming efficiency

The power expended by the swimmer is given by

(4.2)$$\begin{eqnarray}P=\int _{S(t)}\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D748}\boldsymbol{\cdot }\boldsymbol{u}\,\text{d}S=\int _{V(t)}(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748})\boldsymbol{\cdot }\boldsymbol{u}\,\text{d}V+\int _{V(t)}2\unicode[STIX]{x1D707}\unicode[STIX]{x1D640}\boldsymbol{ : }\unicode[STIX]{x1D640}\,\text{d}V.\end{eqnarray}$$

Here, $\unicode[STIX]{x1D748}$ and $\unicode[STIX]{x1D640}$ are the stress tensor and the rate-of-strain tensor, respectively; $S(t)$ represents the surface of the swimmer; and $V$ represents the fluid volume which is bounded by the sheet surface and extends to infinity. The rate-of-strain tensor can be written as $\unicode[STIX]{x1D640}=\frac{1}{2}(\unicode[STIX]{x1D735}\boldsymbol{u}+\unicode[STIX]{x1D735}\boldsymbol{u}^{\text{T}})$, where $\boldsymbol{u}=(u_{1},v_{1})$ represents the first-order velocity field. As inertia is assumed to be negligible in our work, $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}=0$, and the power dissipation is given by

$$\begin{eqnarray}P=\int _{V(t)}2(\unicode[STIX]{x1D707}_{\infty }-\unicode[STIX]{x1D6FE}y+\unicode[STIX]{x1D707}^{\prime })\unicode[STIX]{x1D640}\boldsymbol{ : }\unicode[STIX]{x1D640}\,\text{d}V.\end{eqnarray}$$

Using the expression for the first-order velocity and viscosity fields derived in (3.2) and (3.3), we find that

$$\begin{eqnarray}\int _{V(t)}\unicode[STIX]{x1D707}_{1}^{\prime }\unicode[STIX]{x1D640}\boldsymbol{ : }\unicode[STIX]{x1D640}\,\text{d}V=0.\end{eqnarray}$$

We can now split the power expenditure into two components, $P_{h}=2\unicode[STIX]{x1D707}_{\infty }\int _{V(t)}\unicode[STIX]{x1D640}\boldsymbol{ : }\unicode[STIX]{x1D640}\,\text{d}V$ and $P_{v}=2\unicode[STIX]{x1D6FE}\int _{V(t)}y\unicode[STIX]{x1D640}\boldsymbol{ : }\unicode[STIX]{x1D640}\,\text{d}V$ which denote the magnitude of the mean power expended by the swimmer in a homogeneous fluid of viscosity $\unicode[STIX]{x1D707}_{\infty }$ and the contribution of the ambient viscosity gradient to the power expenditure, respectively. We substitute the expression for the first-order velocity fields and evaluate the integrals to obtain the following expression:

(4.3)$$\begin{eqnarray}P_{v}=\unicode[STIX]{x1D6FC}\unicode[STIX]{x03C0}P_{h}.\end{eqnarray}$$

We find that the viscosity-gradient-induced power expenditure increases with $\unicode[STIX]{x1D6FC}$. This is expected, as the swimmer dissipates more energy in the fluid for larger ambient viscosity gradients.

Lighthill (Reference Lighthill1952) defined swimming efficiency as the ratio of external power required to move the swimmer at a constant velocity to the power expended by the swimmer. This definition is commonly used in the literature to represent the hydrodynamic efficiency of undulating sheets (Li & Ardekani Reference Li and Ardekani2015; Dandekar et al. Reference Dandekar, Shaik and Ardekani2019). We use this definition to quantify the hydrodynamic efficiency of the swimmer in a viscosity-stratified fluid. Hence, we obtain

(4.4)$$\begin{eqnarray}\unicode[STIX]{x1D702}\propto \frac{V^{2}}{P}.\end{eqnarray}$$

Here, $V$ is the swimming velocity and $P$ is the mean power expended by the swimmer over one time period. We use this formula to compare the hydrodynamic efficiency of the swimmer in a homogeneous and a viscosity-stratified fluid. Using (3.11) and (4.3), we obtain

(4.5)$$\begin{eqnarray}\frac{\unicode[STIX]{x1D702}}{\unicode[STIX]{x1D702}_{h}}=\frac{4(\frac{1}{2}+\unicode[STIX]{x1D6FC}f(Pe))^{2}}{1+\unicode[STIX]{x1D6FC}\unicode[STIX]{x03C0}}.\end{eqnarray}$$

As $f(Pe)\sim O(0.01)$ (refer to figure 2), the above expression simplifies to

$$\begin{eqnarray}\frac{\unicode[STIX]{x1D702}}{\unicode[STIX]{x1D702}_{h}}\sim \frac{1}{1+\unicode[STIX]{x1D6FC}\unicode[STIX]{x03C0}}.\end{eqnarray}$$

We thus infer that increasing viscosity gradient decreases the hydrodynamic efficiency of the swimmer, as the swimmer expends power at a faster rate compared to a change in the swimming velocity.