2.1 Orientation probability density

The equation governing the probability density $\unicode[STIX]{x1D6FA}(\boldsymbol{p},t)$ for the swimmer orientation, is given by Subramanian & Koch (Reference Subramanian and Koch2009), Nambiar et al. (Reference Nambiar, Nott and Subramanian2017):

(2.1) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}}{\unicode[STIX]{x2202}t}+\frac{1}{\unicode[STIX]{x1D70F}}\left(\unicode[STIX]{x1D6FA}-\int K(\boldsymbol{p}|\boldsymbol{p}^{\prime })\unicode[STIX]{x1D6FA}(\boldsymbol{p}^{\prime },t)\,\text{d}\boldsymbol{p}^{\prime }\right)-D_{r}\unicode[STIX]{x1D6FB}_{\boldsymbol{p}}^{2}\unicode[STIX]{x1D6FA}=-\unicode[STIX]{x1D735}_{\boldsymbol{ p}}\boldsymbol{\cdot }(\dot{\boldsymbol{p}}\unicode[STIX]{x1D6FA}). & & \displaystyle\end{eqnarray}$$

The terms within the brackets in (2.1) capture the evolution of the orientation distribution due to run-and-tumble dynamics; $K(\boldsymbol{p}|\boldsymbol{p}^{\prime })=\unicode[STIX]{x1D6FD}/(4\unicode[STIX]{x03C0}\sinh \unicode[STIX]{x1D6FD})\exp (\boldsymbol{p}\boldsymbol{\cdot }\boldsymbol{p}^{\prime })$ denotes the correlation between the pre- and post-tumble orientations, with $\unicode[STIX]{x1D6FD}$ being the correlation parameter that is zero for random tumbling, and is approximately 1 for typical wild-type E. coli motion on the micro-scale (Berg Reference Berg1993; Subramanian & Koch Reference Subramanian and Koch2009). The orientational Laplacian arises from the (athermal) rotary diffusion of the bacterium with diffusivity $D_{r}$ . The effect of rotation due to an ambient linear flow is captured by the term on the right side of (2.1), where $\dot{\boldsymbol{p}}$ for slender bodies in an ambient linear flow is given by Jeffery (Reference Jeffery1922)

(2.2) $$\begin{eqnarray}\displaystyle \dot{\boldsymbol{p}}=\unicode[STIX]{x1D74E}\boldsymbol{\cdot }\boldsymbol{p}+\boldsymbol{E}\boldsymbol{\cdot }\boldsymbol{p}-(\boldsymbol{E}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p})\boldsymbol{p}. & & \displaystyle\end{eqnarray}$$

Here, $\unicode[STIX]{x1D74E}$ and $\boldsymbol{E}$ are the vorticity and rate of strain tensors of the ambient flow field, respectively. In (2.1) the spatial dependence of $\unicode[STIX]{x1D6FA}$ has been neglected, since we consider imposed homogeneous linear flows; the bacterial concentration is assumed to be small enough for the effects of the bacterial stress on the imposed flow to be neglected. Interestingly, accounting for effects of the bacterial stress on the imposed flow leads to a homogeneous sheared suspension of pushers becoming linearly unstable to layering perturbations (Laxminarsimharao, Garg & Subramanian Reference Laxminarsimharao, Garg and Subramanian2018). A more detailed description of the various terms and their significance is given in Nambiar et al. (Reference Nambiar, Nott and Subramanian2017).

We cast (2.1) in dimensionless form by scaling time with $\unicode[STIX]{x1D70F}$ , and $\boldsymbol{E}$ , $\unicode[STIX]{x1D74E}$ with $\dot{\unicode[STIX]{x1D6FE}_{0}}$ (the characteristic strength of the imposed velocity gradient), but retain the original notation for brevity, to get

(2.3) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}}{\unicode[STIX]{x2202}t}+\left(\unicode[STIX]{x1D6FA}-\int K(\boldsymbol{p}|\boldsymbol{p}^{\prime })\unicode[STIX]{x1D6FA}(\boldsymbol{p}^{\prime },t)\,\text{d}\boldsymbol{p}^{\prime }\right)-(D_{r}\unicode[STIX]{x1D70F})\unicode[STIX]{x1D6FB}_{\boldsymbol{p}}^{2}\unicode[STIX]{x1D6FA}=-Pe\unicode[STIX]{x1D735}_{\boldsymbol{ p}}\boldsymbol{\cdot }(\dot{\boldsymbol{p}}\unicode[STIX]{x1D6FA}). & & \displaystyle\end{eqnarray}$$

Nambiar et al. (Reference Nambiar, Nott and Subramanian2017) solved (2.3) for $Pe\ll 1$ , as a regular perturbation expansion in $Pe$ keeping only the leading term to characterize the linear response. In appendix A, we extend their solution to higher orders in $Pe$ , and obtain the probability density to $O(Pe^{2})$ for axisymmetric extensional flow, and to $O(Pe^{3})$ for simple shear. In simple shear, the higher-order terms of even powers of $Pe$ contribute corrections to the normal stress differences, whereas, the terms of odd power add corrections to the shear viscosity. Thus, an analysis to $O(Pe^{3})$ is required to capture the leading-order shear-rate dependence of the viscosity (Brenner Reference Brenner1974). In uniaxial extension, however, all the higher-order terms only add corrections to the extensional viscosity, and therefore, an $O(Pe^{2})$ correction suffices. For arbitrary $Pe$ , we resort to a numerical solution using the Galerkin method (Scheraga, Edsall & Gadd Reference Scheraga, Edsall and Gadd1951; Scheraga Reference Scheraga1955; Doi & Edwards Reference Doi and Edwards1978; Strand, Kim & Karrila Reference Strand, Kim and Karrila1987; Chen & Koch Reference Chen and Koch1996). The method is based on expressing $\unicode[STIX]{x1D6FA}$ as a series in spherical harmonics

(2.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FA}(\boldsymbol{p},t)=\mathop{\sum }_{n=0}^{\infty }\mathop{\sum }_{m=-n}^{n}a_{nm}(t)Y_{n}^{m}(\boldsymbol{p}), & & \displaystyle\end{eqnarray}$$

where, $Y_{n}^{m}$ are the surface spherical harmonics defined by

(2.5) $$\begin{eqnarray}\displaystyle Y_{n}^{m}(\boldsymbol{p})\equiv Y_{n}^{m}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})=\sqrt{(2n+1)/(4\unicode[STIX]{x03C0})(n-m)!/(n+m)!}P_{n}^{m}(\cos \unicode[STIX]{x1D703})\exp (\text{i}m\unicode[STIX]{x1D719}) & & \displaystyle\end{eqnarray}$$

and $P_{n}^{m}$ are the associated Legendre Polynomials (Abramowitz & Stegun Reference Abramowitz and Stegun1970). The integral kernel in (2.3) can then be expressed in terms of the Legendre polynomials as

(2.6) $$\begin{eqnarray}\displaystyle K(\boldsymbol{p}|\boldsymbol{p}^{\prime })=\mathop{\sum }_{n=0}^{\infty }A_{n}P_{n}(\boldsymbol{p}\boldsymbol{\cdot }\boldsymbol{p}^{\prime })=\mathop{\sum }_{n=0}^{\infty }\mathop{\sum }_{m,=-n}^{n}\frac{4\unicode[STIX]{x03C0}A_{n}}{(2n+1)}Y_{n}^{m}(\boldsymbol{p})Y_{n}^{m\ast }(\boldsymbol{p}^{\prime }), & & \displaystyle\end{eqnarray}$$

where we have used the addition theorem of spherical harmonics, and the coefficients $A_{n}$ are obtained from the orthogonality relations given in (A 5). Note that for random-tumbling swimmers $A_{n}=0~\forall n\geqslant 1$ .

For axisymmetric extension, $m=0$ , and the expansion reduces to one in Legendre polynomials, $\unicode[STIX]{x1D6FA}(\boldsymbol{p},t)=\sum _{n}a_{n}(t)P_{n}(\cos \unicode[STIX]{x1D703})$ , with $\unicode[STIX]{x1D703}$ being the polar angle measured from the extensional axis. For simple shear flow, $\unicode[STIX]{x1D703}$ and $\unicode[STIX]{x1D719}$ correspond to the polar angle measured from the gradient axis towards the flow–vorticity plane and the azimuthal angle measured from the flow axis in the flow–vorticity plane, respectively (Doi & Edwards Reference Doi and Edwards1978). It will be seen later that this rather unconventional choice of the polar axis for simple shear allows one to express the flow-induced rotation term in (2.3) in a finite number of terms involving the Wigner 3j symbols. In both expansions above, the leading constants ( $a_{0}$ for uniaxial extension, and $a_{00}$ for simple shear) are directly obtained from the normalization constraint for the orientation probability density: $\int \unicode[STIX]{x1D6FA}(\boldsymbol{p},t)\,\text{d}\boldsymbol{p}=1$ . The numerical solution is further complemented by a large- $Pe$ analysis in § 4. We now proceed to solve (2.3) for axisymmetric extensional and simple shear flows.

2.1.1 Axisymmetric extensional flow

The velocity field for axisymmetric extensional flow is $\boldsymbol{u}=(-x_{1}/2,-x_{2}/2,x_{3})$ . Equation (2.2) then allows one to calculate the rate of rotation of the bacterium based on the imposed extensional flow. Substituting the expansion for the orientation distribution, the rate of rotation by the imposed flow and the expansion for $K(\boldsymbol{p}|\boldsymbol{p}^{\prime })$ in (2.3), we have:

(2.7) $$\begin{eqnarray}\displaystyle & & \displaystyle \mathop{\sum }_{n=0}^{\infty }\left[\frac{\text{d}a_{n}}{\text{d}t}+a_{n}-\frac{4\unicode[STIX]{x03C0}A_{n}}{(2n+1)}a_{n}+(D_{r}\unicode[STIX]{x1D70F})n(n+1)a_{n}\right]P_{n}\nonumber\\ \displaystyle & & \displaystyle \quad =\frac{3}{2}\frac{Pe\dot{\unicode[STIX]{x1D6FE}}}{\sin \unicode[STIX]{x1D703}}\frac{\text{d}}{\text{d}\unicode[STIX]{x1D703}}\left(\sin ^{2}\unicode[STIX]{x1D703}\cos \unicode[STIX]{x1D703}\mathop{\sum }_{n=0}^{\infty }a_{n}P_{n}\right).\end{eqnarray}$$

Here, $\dot{\unicode[STIX]{x1D6FE}}~(\equiv \dot{\unicode[STIX]{x1D6FE}}(t))$ is the non-dimensional time dependent rate of extension. Expressing the trigonometric functions in terms of the Legendre polynomials, and using the latter’s recursion and orthogonality relations, results in the following system of linear ordinary differential equations for the series coefficients:

(2.8) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\text{d}a_{n}}{\text{d}t}+a_{n}-\frac{4\unicode[STIX]{x03C0}A_{n}}{(2n+1)}a_{n}+(D_{r}\unicode[STIX]{x1D70F})n(n+1)a_{n}=(2n+1)Pe\dot{\unicode[STIX]{x1D6FE}}\nonumber\\ \displaystyle & & \displaystyle \quad \times \,\mathop{\sum }_{n^{\prime }=0}^{\infty }a_{n^{\prime }}\left[\frac{3}{2}\langle n|2|n^{\prime }\rangle +\frac{3}{4}(n^{\prime }+1)\langle n|1|n^{\prime }+1\rangle -\frac{(n^{\prime }+1)}{4}(\langle n|0|n^{\prime }\rangle +2\langle n|2|n^{\prime }\rangle )\right],\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where the integral of the product of three Legendre polynomials is written compactly in terms of the Wigner 3j symbols (Messiah Reference Messiah1958; Arfken & Weber Reference Arfken and Weber2013),

(2.9) $$\begin{eqnarray}\displaystyle \int _{-1}^{1}P_{l}(x)P_{m}(x)P_{n}(x)\,\text{d}x\equiv \langle l|m|n\rangle =2\left(\begin{array}{@{}ccc@{}}l & m & n\\ 0 & 0 & 0\end{array}\right)\left(\begin{array}{@{}ccc@{}}l & m & n\\ 0 & 0 & 0\end{array}\right). & & \displaystyle\end{eqnarray}$$

2.1.2 Simple shear flow

The flow field for simple shear is $\boldsymbol{u}=(x_{2},0,0)$ , whence the only non-vanishing elements of $\boldsymbol{E}$ and $\unicode[STIX]{x1D714}$ are $E_{12}=E_{21}=\dot{\unicode[STIX]{x1D6FE}}/2$ and $\unicode[STIX]{x1D714}_{12}=-\unicode[STIX]{x1D714}_{21}=\dot{\unicode[STIX]{x1D6FE}}/2$ , $\dot{\unicode[STIX]{x1D6FE}}\equiv \dot{\unicode[STIX]{x1D6FE}}(t)$ now being the non-dimensional time dependent shear rate. Following the same approach as for axisymmetric extension, the system of equations for the $a_{lm}$ is:

(2.10) $$\begin{eqnarray}\displaystyle & & \displaystyle \mathop{\sum }_{l,m}\left[\frac{\text{d}a_{lm}}{\text{d}t}+\left(a_{lm}-a_{lm}\frac{4\unicode[STIX]{x03C0}A_{l}}{(2l+1)}\right)-(D_{r}\unicode[STIX]{x1D70F})l(l+1)a_{lm}\right]Y_{l}^{m}\nonumber\\ \displaystyle & & \displaystyle \quad =-Pe\dot{\unicode[STIX]{x1D6FE}}\left[\left(\frac{2}{3}\sqrt{\frac{4\unicode[STIX]{x03C0}}{5}}Y_{2}^{0}+\frac{\sqrt{4\unicode[STIX]{x03C0}}}{3}Y_{0}^{0}\right)\mathop{\sum }_{l,m}a_{lm}\text{i}L_{y}|Y_{l}^{m}\rangle \right.\nonumber\\ \displaystyle & & \displaystyle \left.\qquad +\,\sqrt{\frac{2\unicode[STIX]{x03C0}}{15}}(Y_{2}^{-1}+Y_{2}^{1})\mathop{\sum }_{l,m}a_{lm}L_{z}|Y_{l}^{m}\rangle \right.\nonumber\\ \displaystyle & & \displaystyle \left.\qquad -\,3\sqrt{\frac{2\unicode[STIX]{x03C0}}{15}}(Y_{2}^{-1}-Y_{2}^{1})\mathop{\sum }_{l,m}a_{lm}Y_{l}^{m}\right].\end{eqnarray}$$

The definitions of the angular momentum operators in the above equation ( $L_{y}\{\cdot \}$ and $L_{z}\{\cdot \}$ ) are given by (see appendix of Doi & Edwards (Reference Doi and Edwards1978)) (Messiah Reference Messiah1958):

(2.11) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}rcl@{}}\text{i}L_{y}|l,m\rangle & = & \frac{1}{2}\sqrt{l(l+1)-m(m+1)}|l,m+1\rangle \\ & & -\,\frac{1}{2}\sqrt{l(l+1)-m(m-1)}|l,m-1\rangle \\ L_{z}|l,m\rangle & = & m|l,m\rangle .\end{array}\right\} & & \displaystyle\end{eqnarray}$$

Finally, using the orthogonality of the spherical harmonics, one obtains:

(2.12) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\text{d}a_{lm}}{\text{d}t}+a_{lm}-a_{lm}\frac{4\unicode[STIX]{x03C0}A_{l}}{(2l+1)}-(D_{r}\unicode[STIX]{x1D70F})l(l+1)a_{lm}=-Pe\dot{\unicode[STIX]{x1D6FE}}\nonumber\\ \displaystyle & & \displaystyle \quad \times \,\mathop{\sum }_{l^{\prime },m^{\prime }}\left[\left\{\frac{1}{3}\sqrt{\frac{4\unicode[STIX]{x03C0}}{5}}(a_{l^{\prime }m^{\prime }}^{\prime }\langle l,m|Y_{2}^{0}|l^{\prime },m^{\prime }+1\rangle -a_{l^{\prime }m^{\prime }}^{\prime \prime }\langle l,m|Y_{2}^{0}|l^{\prime },m^{\prime }-1\rangle )\right.\right.\nonumber\\ \displaystyle & & \displaystyle \quad +\left.\frac{\sqrt{4\unicode[STIX]{x03C0}}}{6}(a_{l^{\prime }m^{\prime }}^{\prime }\langle l,m|Y_{0}^{0}|l^{\prime },m^{\prime }+1\rangle +a_{l^{\prime }m^{\prime }}^{\prime \prime }\langle l,m|Y_{0}^{0}|l^{\prime },m^{\prime }-1\rangle )\right\}\nonumber\\ \displaystyle & & \displaystyle \quad +\left.\sqrt{\frac{2\unicode[STIX]{x03C0}}{15}}\{a_{l^{\prime }m^{\prime }}(\langle l,m|Y_{2}^{-1}|l^{\prime },m^{\prime }\rangle (m^{\prime }+3)+\langle l,m|Y_{2}^{1}|l^{\prime },m^{\prime }\rangle (m^{\prime }-3))\}\right],\end{eqnarray}$$

where, $a_{l^{\prime }m^{\prime }}^{\prime }=a_{l^{\prime }m^{\prime }}\sqrt{l^{\prime }(l^{\prime }+1)-m^{\prime }(m^{\prime }+1)}$ , $a_{l^{\prime }m^{\prime }}^{\prime \prime }=a_{l^{\prime }m^{\prime }}\sqrt{l^{\prime }(l^{\prime }+1)-m^{\prime }(m^{\prime }-1)}$ and $\langle l,m|Y_{p}^{q}|l^{\prime },m^{\prime }\rangle$ is given by Doi & Edwards (Reference Doi and Edwards1978):

(2.13) $$\begin{eqnarray}\displaystyle & & \displaystyle \langle l,m|Y_{p}^{q}|l^{\prime },m^{\prime }\rangle \nonumber\\ \displaystyle & & \displaystyle \quad =(-1)^{m}\left[\frac{(2l^{\prime }+1)(2l+1)(2p+1)}{4\unicode[STIX]{x03C0}}\right]^{1/2}\left(\begin{array}{@{}ccc@{}}l & p & l^{\prime }\\ -m & q & m^{\prime }\end{array}\right)\left(\begin{array}{@{}ccc@{}}l & p & l^{\prime }\\ 0 & 0 & 0\end{array}\right).\end{eqnarray}$$

We solve the systems (2.8) and (2.12) numerically using the biorthogonal expansion technique (Strand et al. Reference Strand, Kim and Karrila1987). For simple shear flow at $Pe=100$ , $l_{Max}=64$ corresponding to a total of 4225 terms in the spherical harmonics series. For axisymmetric extension, $l_{Max}=264$ , which is also the actual number of terms.

2.2 Bulk hydrodynamic stress

The deviatoric part of the bacterial contribution to the bulk hydrodynamic stress can be obtained as an orientation space average of the stresslet $\boldsymbol{S}(\boldsymbol{p})$ as:

(2.14) $$\begin{eqnarray}\displaystyle \langle \unicode[STIX]{x1D748}^{b}\rangle =n\int \text{d}\boldsymbol{p}\unicode[STIX]{x1D6FA}(\boldsymbol{p},t)\boldsymbol{S}(\boldsymbol{p}), & & \displaystyle\end{eqnarray}$$

where $n$ is the bacterial number density. Using the definition of the stresslet in terms of the swimmer force density, we obtain

(2.15) $$\begin{eqnarray}\displaystyle \langle \unicode[STIX]{x1D748}^{b}\rangle =n\int \text{d}\boldsymbol{p}\unicode[STIX]{x1D6FA}(\boldsymbol{p},t)\int _{-L/2}^{L/2}s\left[\boldsymbol{p}\boldsymbol{f}(s)+\boldsymbol{f}(s)\boldsymbol{p}-\frac{2}{3}(\boldsymbol{f}(s)\boldsymbol{\cdot }\boldsymbol{p})\boldsymbol{I}\right]\,\text{d}s, & & \displaystyle\end{eqnarray}$$

where the inner integral is over the coordinate $s$ along the axis of the slender swimmer having total length (head $+$ tail) $L$ , and $\boldsymbol{f}(s)$ is the linear force density. The force density has two contributions, one arising due to the inextensibility of the swimming bacterium (passive) and another due to its intrinsic activity, that is, $\boldsymbol{f}(s)=\boldsymbol{f}^{p}(s)+\boldsymbol{f}^{a}(s)$ . Substituting the expressions for the two contributions, as given in Nambiar et al. (Reference Nambiar, Nott and Subramanian2017), the deviatoric part of the bacterial stress scaled by $\unicode[STIX]{x1D707}\dot{\unicode[STIX]{x1D6FE}}$ reduces to:

(2.16) $$\begin{eqnarray}\displaystyle \langle \unicode[STIX]{x1D748}^{b}\rangle & = & \displaystyle \frac{\unicode[STIX]{x03C0}(nL^{3})}{6\ln \unicode[STIX]{x1D705}}\dot{\unicode[STIX]{x1D6FE}}(t)\int \text{d}\boldsymbol{p}\,\unicode[STIX]{x1D6FA}(\boldsymbol{p},t)(\boldsymbol{E}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p})\left(\boldsymbol{p}\boldsymbol{p}-\frac{\boldsymbol{I}}{3}\right)\nonumber\\ \displaystyle & & \displaystyle -\,\frac{\unicode[STIX]{x1D6FC}}{2M}\frac{nUL^{2}\unicode[STIX]{x1D70F}}{Pe}\int \text{d}\boldsymbol{p}\unicode[STIX]{x1D6FA}(\boldsymbol{p},t)\left(\boldsymbol{p}\boldsymbol{p}-\frac{\boldsymbol{I}}{3}\right),\end{eqnarray}$$

where $\unicode[STIX]{x1D705}$ is the aspect ratio of the slender swimmer ( $\unicode[STIX]{x1D705}\gg 1$ ). The first term on the right side of the above equation is the passive hydrodynamic contribution, and the second term is the active hydrodynamic contribution, due to the intrinsic force dipoles, with $M$ being the head mobility coefficient and $\unicode[STIX]{x1D6FC}$ being the ratio of the head to total bacterium length.

2.2.1 Axisymmetric extensional flow

On account of axisymmetry, determination of one of the diagonal stress components, $\langle \unicode[STIX]{x1D70E}_{E}^{b}\rangle _{33}$ from (2.16), is sufficient, and this takes the form:

(2.17) $$\begin{eqnarray}\displaystyle \langle \unicode[STIX]{x1D70E}_{E}^{b}\rangle _{33} & = & \displaystyle \frac{\unicode[STIX]{x03C0}(nL^{3})}{9\ln \unicode[STIX]{x1D705}}\dot{\unicode[STIX]{x1D6FE}}(t)\int \text{d}\boldsymbol{p}\unicode[STIX]{x1D6FA}(\boldsymbol{p},t)(P_{2}(\cos \unicode[STIX]{x1D703}))^{2}\nonumber\\ \displaystyle & & \displaystyle -\,\frac{\unicode[STIX]{x1D6FC}}{3M}\frac{nUL^{2}\unicode[STIX]{x1D70F}}{Pe}\int \text{d}\boldsymbol{p}\unicode[STIX]{x1D6FA}(\boldsymbol{p},t)P_{2}(\cos \unicode[STIX]{x1D703}).\end{eqnarray}$$

Using the orthogonality of the Legendre polynomials, the orientation integrals in (2.17) may be written in terms of the coefficients $a_{n}$ , the result being:

(2.18) $$\begin{eqnarray}\displaystyle \langle \unicode[STIX]{x1D70E}_{E}^{b}\rangle _{33}=\frac{4\unicode[STIX]{x03C0}^{2}(nL^{3})}{45\ln \unicode[STIX]{x1D705}}\dot{\unicode[STIX]{x1D6FE}}(t)\left[a_{0}(t)+\frac{2}{7}(a_{2}(t)+a_{4}(t))\right]-\frac{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FC}(nUL^{2}\unicode[STIX]{x1D70F})}{15MPe}a_{2}(t), & & \displaystyle\end{eqnarray}$$

with $a_{n}(t)$ being the coefficients defined in § 2.1.1.

2.2.2 Simple shear flow

Here, the lone off-diagonal term, viz. $\langle \unicode[STIX]{x1D70E}_{S}^{b}\rangle _{12}$ , gives the shear stress, while the diagonal terms yield the normal stress differences. Proceeding in a manner similar to (2.2.1), the shear viscosity is expressed as

(2.19) $$\begin{eqnarray}\displaystyle \langle \unicode[STIX]{x1D70E}_{S}^{b}\rangle _{12} & = & \displaystyle \frac{\unicode[STIX]{x03C0}(nL^{3})}{6\ln \unicode[STIX]{x1D705}}\dot{\unicode[STIX]{x1D6FE}}(t)\left[\frac{2\sqrt{\unicode[STIX]{x03C0}}}{15}a_{0,0}(t)+\frac{2}{21}\sqrt{\frac{\unicode[STIX]{x03C0}}{5}}a_{2,0}(t)+\frac{1}{7}\sqrt{\frac{2\unicode[STIX]{x03C0}}{15}}(a_{2,2}(t)+a_{2,-2}(t))\right.\nonumber\\ \displaystyle & & \displaystyle -\left.\frac{8\sqrt{\unicode[STIX]{x03C0}}}{105}a_{4,0}(t)+\frac{2}{21}\sqrt{\frac{2\unicode[STIX]{x03C0}}{5}}(a_{4,2}(t)+a_{4,-2}(t))\right]\nonumber\\ \displaystyle & & \displaystyle +\,\frac{\unicode[STIX]{x1D6FC}(nUL^{2}\unicode[STIX]{x1D70F})}{2MPe}\sqrt{\frac{2\unicode[STIX]{x03C0}}{15}}\left[a_{2,1}(t)-a_{2,-1}(t)\right].\end{eqnarray}$$

Similarly, the normal stress differences, $N_{1}=(\langle \unicode[STIX]{x1D70E}_{S}^{b}\rangle _{11}-\langle \unicode[STIX]{x1D70E}_{S}^{b}\rangle _{22})/(\dot{\unicode[STIX]{x1D6FE}}\unicode[STIX]{x1D70F})$ ( $\equiv (\unicode[STIX]{x1D70E}_{xx}-\unicode[STIX]{x1D70E}_{yy})/(\unicode[STIX]{x1D707}\dot{\unicode[STIX]{x1D6FE}}^{2}\unicode[STIX]{x1D70F})$ ) and $N_{2}=(\langle \unicode[STIX]{x1D70E}_{S}^{b}\rangle _{22}-\langle \unicode[STIX]{x1D70E}_{S}^{b}\rangle _{33})/(\dot{\unicode[STIX]{x1D6FE}}\unicode[STIX]{x1D70F})$ ( $\equiv (\unicode[STIX]{x1D70E}_{yy}-\unicode[STIX]{x1D70E}_{zz})/(\unicode[STIX]{x1D707}\dot{\unicode[STIX]{x1D6FE}}^{2}\unicode[STIX]{x1D70F})$ ) are defined as

(2.20a ) $$\begin{eqnarray}\displaystyle N_{1} & = & \displaystyle \frac{\unicode[STIX]{x03C0}(nL^{3})}{6\ln \unicode[STIX]{x1D705}}\frac{1}{Pe}\dot{\unicode[STIX]{x1D6FE}}(t)\left[\frac{1}{3}\sqrt{\frac{\unicode[STIX]{x03C0}}{5}}(a_{4,1}(t)-a_{4,-1}(t))-\frac{1}{3}\sqrt{\frac{\unicode[STIX]{x03C0}}{35}}(a_{4,3}(t)-a_{4,-3}(t))\right]\nonumber\\ \displaystyle & & \displaystyle +\,\frac{\unicode[STIX]{x1D6FC}(nUL^{2}\unicode[STIX]{x1D70F})}{2MPe^{2}}\left[2\sqrt{\frac{\unicode[STIX]{x03C0}}{5}}a_{2,0}(t)-\sqrt{\frac{2\unicode[STIX]{x03C0}}{15}}(a_{2,2}(t)+a_{2,-2}(t))\right],\end{eqnarray}$$

(2.20b ) $$\begin{eqnarray}\displaystyle N_{2} & = & \displaystyle \frac{\unicode[STIX]{x03C0}(nL^{3})}{6\ln \unicode[STIX]{x1D705}}\frac{1}{Pe}\dot{\unicode[STIX]{x1D6FE}}(t)\left[\frac{2}{7}\sqrt{\frac{2\unicode[STIX]{x03C0}}{15}}(a_{2,-1}(t)-a_{2,1}(t))+\frac{5}{21}\sqrt{\frac{\unicode[STIX]{x03C0}}{5}}(a_{4,-1}(t)-a_{4,1}(t))\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\frac{1}{3}\sqrt{\frac{\unicode[STIX]{x03C0}}{35}}(a_{4,-3}(t)-a_{4,3}(t))\right]\nonumber\\ \displaystyle & & \displaystyle -\,\frac{\unicode[STIX]{x1D6FC}(nUL^{2}\unicode[STIX]{x1D70F})}{2MPe^{2}}\left[2\sqrt{\frac{\unicode[STIX]{x03C0}}{5}}a_{2,0}(t)+\sqrt{\frac{2\unicode[STIX]{x03C0}}{15}}(a_{2,2}(t)+a_{2,-2}(t))\right].\end{eqnarray}$$

The expressions derived in §§ 2.2.1 and 2.2.2 are valid for an arbitrary shear-rate history, but in what follows, we restrict ourselves to a study of the response in step shear.

3.1 Step initiation of the imposed flow

Here, we examine the response of the suspension to impulsive initiation of the flow, implying $\dot{\unicode[STIX]{x1D6FE}}(t)=H(t)$ in (2.18)–(2.20). For axisymmetric extension, at any $Pe$ , the flow distorts the isotropic distribution, creating an excess of orientations in the vicinity of the extensional axis. Figure 1(a,b) depicts the iso-contours of the orientation probability in extensional flow at the initial instant, an intermediate time and at steady state, for $Pe=0.1$ and 10, respectively. In figure 2, we plot the bacterial contribution to the time dependent extensional viscosity for different $Pe$ . For $Pe\ll 1$ , the approach towards the anisotropic steady state occurs in a time of $O(\unicode[STIX]{x1D70F})$ , and the oriented extensile swimmer dipoles reinforce the imposed extension, implying a reduction in the total viscosity (below that of the solvent). However, with increasing $Pe$ , the passive contribution to the stress begins to increase and is eventually dominant, implying a steady-state extensional viscosity higher than that of the solvent. Therefore, the steady-state extensional viscosity plateau, starts off being negative at small $Pe$ , passes through a minimum at $Pe\sim 0.3$ and changes sign at $Pe\sim 10$ . Although the said plateau is positive at $Pe=10$ , implying a dominant passive stress component, the relaxation to this steady state is nevertheless non-monotonic, implying a non-trivial effect of activity on the transient even at large $Pe$ (see inset of figure 2). For $Pe=50$ or larger, however, the active component has a negligible effect at all times, and a monotonic relaxation prevails. Note that as $Pe$ increases from 0.01 to 50, in figure 2, the characteristic relaxation time reduces from being $O(\unicode[STIX]{x1D70F})$ when $Pe\ll 1$ to $O(\dot{\unicode[STIX]{x1D6FE}}^{-1})$ when $Pe\gg 1$ , as expected.

Figure 2.

The time dependent bacterial extensional viscosity ( $\langle \unicode[STIX]{x1D70E}_{E}^{b}\rangle _{33}$ ), given by (2.18) for step initiation of the flow at $t=0$ , plotted for $Pe\sim 0.01{-}50$ . The inset consists of the extensional viscosity for $Pe=10$ and $Pe=50$ , with the time being scaled by the inverse rate of extension.

Figure 3(a,b) shows the iso-contours of the orientation probability in simple shear, again at the initial instant, an intermediate time and at steady state, for $Pe=0.022$ and 63.9. At $Pe=0.022$ , there is an enhanced orientation probability along the extensional axis (at $45^{\circ }$ to the flow direction) that dictates the linear response (see below). In contrast, at $Pe=63.9$ , the peak in the probability density is rotated by the shear into near alignment with the flow direction. Figure 4 contains only the bacterial contribution to the time dependent viscosity of the swimmer suspension in response to an imposed step initiation of simple shear (the solvent contribution is unity). The collapse of the relaxation curves for $Pe$ up to 0.075 is indicative of a linear response regime ( $Pe\ll 1$ ). An explanation similar to that for extensional flow above can be adapted to simple shear, where the dipoles associated with the extensional-axis-aligned swimmers lead to a reduced viscosity (Nambiar et al. Reference Nambiar, Nott and Subramanian2017). With increasing $Pe$ , there is a transition to a nonlinear response arising from the rotation of the dipoles away from the extensional axis. This misalignment leads to a shear-thickening behaviour until $Pe\approx 10$ (it is shown in appendix A that a small- $Pe$ expansion of the orientation probability density, at $O(Pe^{3})$ , captures the shear-rate dependence of the viscosity until $Pe\approx 0.3$ ). Although the steady viscosity is shear thickening for $Pe\in [3,10]$ , a non-monotonic transient in the viscosity is observed, that is, the viscosity reaches a minimum before rising to a plateau higher than that in the linear regime. This is due to the active component, which continues to remains important at moderate $Pe$ , as the bacteria are transiently aligned with the extensional axis on (short) time scales of $O(\dot{\unicode[STIX]{x1D6FE}}^{-1})$ . Further increase in $Pe$ leads to a transition, from an active stress-dominated regime with a shear-thickening rheology at low to moderate $Pe$ , to a passive stress-dominated regime, which leads to an eventual shear-thinning behaviour (see § 4.1.2 for more details). The latter regime is similar to that reported by Strand et al. (Reference Strand, Kim and Karrila1987) for a suspension of passive rigid rods. This transition occurs over an interval of $Pe$ between 10 and 30, and a manifestation of this transition is the emergence of a peak for short times owing to the passive hydrodynamic contribution associated with swimmers momentarily aligned with the extensional axis (see inset of figure 4). The characteristic relaxation time reduces from $O(\unicode[STIX]{x1D70F})$ for $Pe\ll 1$ to $O(\dot{\unicode[STIX]{x1D6FE}}^{-1})$ when $Pe\gg 1$ , as in extensional flow. The reader is directed to § 4 for further discussion on the $Pe\gg 1$ regime.

Figure 3.

The iso-contours of the orientation probability density, at three different time instants, plotted for simple shear at (a) $Pe=0.1$ and (b) at $Pe=63.9$ .

Figure 4.

The time dependent shear viscosity ( $\langle \unicode[STIX]{x1D70E}_{S}\rangle _{12}$ ), given by (2.19) for step initiation of the flow at $t=0$ , plotted for $Pe\sim 0.022{-}100$ . The inset consists of the shear viscosity, for $Pe=34.6$ , $Pe=63.9$ and $Pe=100$ , plotted as a function of the time scaled by the inverse shear rate.

3.2 Step cessation of the imposed flow

We now consider the case of impulsive shear cessation, so $\dot{\unicode[STIX]{x1D6FE}}(t)=1-H(t)$ in (2.18)–(2.20). In the absence of an imposed flow, the passive hydrodynamic stress component is absent; the flow merely acts to set up a $Pe$ -dependent initial orientation distribution. Thus, the relaxation of the active stress component remains the same for all $Pe$ , and for both simple shear and axisymmetric extension. This is readily seen from (2.12) by setting the right-hand side to zero. The relaxation dynamics of the orientation distribution is governed by the spectrum corresponding to rotary diffusion and (random) tumbling; thus, the $n$ th spherical harmonic will decay as $\exp [-(n(n+1)\unicode[STIX]{x1D70F}D_{r}+1)t]$ . The active stress being proportional to the second harmonic will decay as $\exp [-(6\unicode[STIX]{x1D70F}D_{r}+1)t]$ regardless of $Pe$ . This is similar to the relaxation response in suspensions of anisotropic Brownian particles (Hinch & Leal Reference Hinch and Leal1973), where the step-down phase only accounts for the Brownian stress that is also proportional to the second harmonic. In figure 5(a) we show the time dependent shear viscosity for the full range of $Pe$ , and relaxations remain qualitatively similar at all $Pe$ . Beginning from an anisotropic distribution of swimmer orientations under steady shear (or extension) at $t=0^{-}$ , run-and-tumble dynamics and rotary diffusion act to randomize the distribution, thereby yielding a zero stress at steady state. Here again the initial conditions following the negative jump are nearly coincident for $Pe$ up to 0.075, indicative of a linear response regime. For $Pe\gg 1$ , the active stress scales as $Pe^{-4/3}$ , and accordingly, the relaxation curves at large $Pe$ collapse on re-scaling the stress with $Pe^{-4/3}$ (see inset of figure 5 a).

For $D_{r}=0.062~\text{s}^{-1}$ and $\unicode[STIX]{x1D70F}=1~\text{s}$ , one obtains a stress relaxation time of $(6\unicode[STIX]{x1D70F}D_{r}+1)$ $\unicode[STIX]{x1D70F}\approx 1.3s$ . This is still less than the observed relaxation times of approximately 2 s. As explained in Nambiar et al. (Reference Nambiar, Nott and Subramanian2017), assuming $\unicode[STIX]{x1D6E4}$ -distributed mean run times, that is, $P_{\unicode[STIX]{x1D6E4}}(\unicode[STIX]{x1D70F};\unicode[STIX]{x1D712},k)=\unicode[STIX]{x1D70F}^{k-1}/(\unicode[STIX]{x1D712}^{k}\unicode[STIX]{x1D6E4}(k))e^{-\unicode[STIX]{x1D70F}/\unicode[STIX]{x1D712}}$ , can account for this underprediction. Here, $k$ is the shape factor and $\unicode[STIX]{x1D712}$ is the scale factor, and the mean $\langle \unicode[STIX]{x1D70F}\rangle \equiv k\unicode[STIX]{x1D712}=1$ is fixed. Figure 5(b) shows the stress relaxations for $Pe=0.022$ and 5, assuming a $\unicode[STIX]{x1D6E4}$ -distributed $\unicode[STIX]{x1D70F}$ with $k=0.62$ , that best fits the observed exponential decay over a time scale of $2\unicode[STIX]{x1D70F}$ .

Figure 5.

The time dependent shear viscosity ( $\langle \unicode[STIX]{x1D70E}_{S}\rangle _{12}$ ), given by (2.19) for step-down from steady flow to rest is plotted (a) for a range of $Pe\sim 0.022{-}100$ . The inset has $\langle \unicode[STIX]{x1D70E}_{S}\rangle _{12}$ scaled by $Pe^{4/3}$ for $Pe=0.138,63.9,100$ . The collapse of the largest two $Pe$ is consistent with the scaling of the active stress in the $Pe\gg 1$ regime. (b) Comparison of the $\unicode[STIX]{x1D6E4}$ distributed relaxation with a shape parameter $k=0.64$ that best fits an exponential relaxation of $2\unicode[STIX]{x1D70F}$ (Nambiar et al. Reference Nambiar, Nott and Subramanian2017) for $Pe=0.022$ and 5.

3.3 Comparison with experiments

Having analysed the stress relaxations in both the step-up and step-down protocols, we now compare our results with the experiments of Lopez et al. (Reference Lopez, Gachelin, Douarche, Auradou and Clement2015). Figure 6(a,b) shows the stress relaxations measured by Lopez et al. (Reference Lopez, Gachelin, Douarche, Auradou and Clement2015) alongside our predictions (with the solvent stress now included) for the same $Pe$ . In accordance with the experimental protocol, the predictions for the step-up and step-down phases are combined, assuming a duration of steady shear for 30 s. Note that this duration is more than sufficient for the initial isotropic microstructure to attain a steady state, and therefore, the combined prediction amounts to a linear superposition of the step-up and step-down phases. There is good agreement for the $Pe$ range corresponding to the linear regime, although the effects of nonlinearity, during the step-up phase, appear at smaller $Pe$ in the experiments. For the shear cessation, the initial stress jump is underpredicted at the highest $Pe$ ( $Pe=63.9$ ). Notwithstanding this disparity, using a $\unicode[STIX]{x1D6E4}$ -distributed $\unicode[STIX]{x1D70F}$ yields stress relaxation times that agree well with the experiments for all $Pe$ . However, there are some important disagreements in the step-up phase that need to be highlighted. While the observations indicate a pronounced minimum in the time dependent viscosity for $Pe=0.471$ and 0.87, this behaviour is observed in our computations only starting at $Pe\approx 3$ (not shown). The more important divergence is the magnitude of the steady-state plateaus. The observations indicate a rapid rise in the plateau values for $Pe>0.075$ (in the nonlinear regime), with the plateau, for $Pe=0.87$ being nearly 82 % of the initial viscous jump. In contrast, the computed relaxation at $Pe=0.87$ , apart from being monotonic, saturates at a value that is only 42 % of the initial viscous jump. More puzzlingly, at the two highest values of $Pe$ (34.6 and 63.9), the observed relaxation is a monotonic rise to a plateau stress that is above the initial viscous jump. This response mimics a suspension of orientable passive Brownian particles at low $Pe$ ; that is, a viscous jump followed by elastic relaxation to a steady state above the viscous Newtonian plateau. This behaviour is contrary to intuition, as at large $Pe$ we expect the flow-aligning effect of the imposed shear (and the resulting decrease in the viscosity) to overwhelm any relaxation to isotropy.

Figure 6.

The time dependent shear viscosity (a) reported in Lopez et al. (Reference Lopez, Gachelin, Douarche, Auradou and Clement2015) (reproduced with permission – © 2015 American Physical Society), and (b) given by (2.19), with the solvent contribution included, plotted for $Pe\sim 0.022{-}63.9$ . For E. coli, the run duration is $\unicode[STIX]{x1D70F}\approx 1~\text{s}$ , and hence $Pe$ equals the shear rate $\dot{\unicode[STIX]{x1D6FE}}$ (in $\text{s}^{-1}$ ).

Another point of divergence between our analysis and the observations of Lopez et al. (Reference Lopez, Gachelin, Douarche, Auradou and Clement2015) is with regard to the relaxation time to the steady-state viscosity during step-up. The relaxation times remain comparable in the experiment for the entire range of $Pe$ examined; it is approximately 3.2 s at low $Pe$ , and slightly above 1 s at the highest $Pe$ . This observation is difficult to explain. In a dilute bacterial suspension (the assumption of diluteness here being reasonable, given that the true volume fraction in Lopez et al. (Reference Lopez, Gachelin, Douarche, Auradou and Clement2015) was around 0.67 %), the primary relaxation is determined by the faster of the two competing effects, namely the stochastic relaxation mechanisms (run-and-tumble dynamics or rotary diffusion) and flow-induced orientation. This implies that, for their largest shear rate of $\dot{\unicode[STIX]{x1D6FE}}=63.9~\text{s}^{-1}$ , the expected relaxation time is 0.016 s. Therefore, the observed relaxation time of approximately 1 s is higher by a factor of 60. Consideration of a distribution of mean run times, did account for slightly larger relaxation times observed in the linear response regime (Nambiar et al. Reference Nambiar, Nott and Subramanian2017), but cannot account for the pronounced discrepancy in the nonlinear response. A possible explanation for the large observed relaxation times is that at the volume fraction considered, the effective volume fraction is $nL^{3}=3.43$ , at which the suspension has transitioned to a regime of collective motion (Krishnamurthy & Subramanian Reference Krishnamurthy and Subramanian2015; Nambiar et al. Reference Nambiar, Nott and Subramanian2017).

3.4 Normal stress differences

Figure 7.

The time dependent bacterial normal stress differences determined numerically in (2.20a ) and (2.20b ) are plotted for $Pe\sim [0.01,100]$ . (a) $N_{1}$ (b) $N_{2}$ .

While the experiments of Lopez et al. (Reference Lopez, Gachelin, Douarche, Auradou and Clement2015) measured only the shear viscosity, we will also look at the normal stress differences for the step-up phase in what follows, since they are representative of the non-trivial elastic signatures associated with pusher suspensions, and might be responsible for triggering elastic instabilities (similar to polymer solutions) in other inhomogeneous flow settings. We first interpret the steady-state results. As shown in appendix B (figure 10 b,c), the first ( $N_{1}$ ) and the second normal stress ( $N_{2}$ ) coefficients begin from a low- $Pe$ plateau (negative for $N_{1}$ and positive for $N_{2}$ ), eventually asymptoting to zero with appropriate power-law scalings for $Pe\rightarrow \infty$ . These trends are consistent with the results of Saintillan (Reference Saintillan2010a ).

The temporal variations of the first and second normal stress differences (NSD) for different $Pe$ are shown in figures 7(a) and 7(b). When $Pe\ll 1$ , the dominant contributions to both $N_{1}$ and $N_{2}$ come from the active component. Therefore, in this limit, the NSD have signs opposite to those expected for passive suspensions ( $N_{1}$ is negative and $N_{2}$ is positive), the relaxation to steady state occurring on a time scale of $O(\unicode[STIX]{x1D70F})$ . With increase in $Pe$ , the time scale for the relaxation of the NSD eventually shrinks to $O(\dot{\unicode[STIX]{x1D6FE}}^{-1})$ . Following the same reasoning as for the shear viscosity, one expects the passive component to be dominant when $Pe\gg 1$ , implying a sign reversal of the NSD for sufficiently large $Pe$ (Brenner Reference Brenner1974; Strand et al. Reference Strand, Kim and Karrila1987). For $Pe=100$ , the peak $N_{1}$ , attained at a time of $O(\dot{\unicode[STIX]{x1D6FE}}^{-1})$ , is indeed positive, despite the eventual plateau being negative. Interestingly, for $N_{2}$ we observe no sign change for the chosen values of the swimmer parameters. Further discussion on the large- $Pe$ scaling behaviour of the NSD is given in appendix B.

3.5 Effect of correlation in tumbling on the viscosity

In (2.3), the tumbling kernel $K(\boldsymbol{p}|\boldsymbol{p}^{\prime })=\unicode[STIX]{x1D6FD}\exp (\unicode[STIX]{x1D6FD}\boldsymbol{p}\boldsymbol{\cdot }\boldsymbol{p}^{\prime })/(4\unicode[STIX]{x03C0}\sinh \unicode[STIX]{x1D6FD})$ (Subramanian & Koch Reference Subramanian and Koch2009), is a function of the parameter $\unicode[STIX]{x1D6FD}$ , going from $1/4\unicode[STIX]{x03C0}$ for $\unicode[STIX]{x1D6FD}=0$ (random tumbling) to $\unicode[STIX]{x1D6FF}(\boldsymbol{p}-\boldsymbol{p}^{\prime })$ for $\unicode[STIX]{x1D6FD}\rightarrow \infty$ . The results shown thus far have all been for random tumbling. To study the role of correlation of adjacent tumbles we consider the cases $\unicode[STIX]{x1D6FD}=0$ and 1, the latter being consistent with the run-and-tumble statistics of wild-type E. coli (Berg Reference Berg1993; Subramanian & Koch Reference Subramanian and Koch2009). The main observation is that correlated tumbling has no effect on the qualitative nature of the stress relaxation (not shown). Quantitatively, correlated tumbling leads to a marginally slower relaxation of the orientation distribution. This, together with the differing dominant contributions to the stress, active at low $Pe$ and passive at large $Pe$ , accounts for the change in the sign of the correlation-induced correction to the random-tumbling plateaus. Thus, correlated tumbles lead to a stronger negative stress at low $Pe$ (until approximately 0.5), and a greater positive stress at higher $Pe$ (greater than approximately 10).

4.1 Outer convective solution ( $Pe\rightarrow \infty$ )

Beginning from an isotropic distribution $\unicode[STIX]{x1D6FA}(\boldsymbol{p},0)=1/(4\unicode[STIX]{x03C0})$ , the evolution of the swimmer orientation probability density is different for axisymmetric extension and simple shear. In the former case, the orientation distribution is stretched towards alignment with the extensional axis. However, in the latter case, the alignment with the extensional axis (at $45^{\circ }$ ) is a transient phenomenon, occurring on time scales of $O(1)$ , following which the orientation distribution is rotated towards the flow axis while also being stretched at the same time. Therefore, in addition to the obvious advantage of azimuthal angle independence (mentioned in § 2), the analysis for axisymmetric extension is also simpler since one only needs to consider stretching of $\unicode[STIX]{x1D6FA}(\boldsymbol{p},t)$ with the peak probability always along the extensional axis. At leading order, one needs to solve for the unsteady convective evolution in (4.1), given by:

(4.2) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{o}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}_{\boldsymbol{p}}\boldsymbol{\cdot }(\dot{\boldsymbol{p}}\unicode[STIX]{x1D6FA}_{o})=0. & & \displaystyle\end{eqnarray}$$

Here, $\unicode[STIX]{x1D6FA}_{o}$ refers to the outer solution.

4.1.1 Axisymmetric extensional flow

For step initiated flow $\dot{\unicode[STIX]{x1D6FE}}=H(t)\equiv 1$ for $t>0$ . Choosing the coordinate system of § 2.1.1, substituting for $\dot{\boldsymbol{p}}$ , and setting $\cos \unicode[STIX]{x1D703}=x$ , simplifies (4.2) to give:

(4.3) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{o}}{\unicode[STIX]{x2202}t}+\frac{3}{2}x(1-x^{2})\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{o}}{\unicode[STIX]{x2202}x}+\frac{3}{2}(1-3x^{2})\unicode[STIX]{x1D6FA}_{o}=0. & & \displaystyle\end{eqnarray}$$

The above equation can be solved by using the method of characteristics. Denoting the factor multiplying $\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{o}/\unicode[STIX]{x2202}x$ as $\text{d}x/\text{d}t$ , one has the following pair of ordinary differential equations (ODE) to solve:

(4.4a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}x}{\text{d}t}=\frac{3}{2}x(1-x^{2}), & \displaystyle\end{eqnarray}$$

(4.4b ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\unicode[STIX]{x1D6FA}_{o}}{\text{d}t}=-\frac{3}{2}(1-3x^{2})\unicode[STIX]{x1D6FA}_{o}, & \displaystyle\end{eqnarray}$$

where $x$ in (4.4b ) is now a function of time. The solution is

(4.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FA}_{o}=\frac{1}{4\unicode[STIX]{x03C0}}\left(\frac{1}{1-\sin ^{2}\unicode[STIX]{x1D703}(1-\exp (3t))}\right)^{3/2}\exp (3t). & & \displaystyle\end{eqnarray}$$

The factor $\exp (3t)$ implies an exponential increase of $\unicode[STIX]{x1D6FA}_{o}(\boldsymbol{p},t)$ along the extensional axis. Conservation of probability implies that the width of the distribution decreases as $\exp (-3t/2)$ , so that, in the absence of stochastic relaxation, $\unicode[STIX]{x1D6FA}_{o}(\boldsymbol{p},t)$ approaches a Dirac delta function for $t\rightarrow \infty$ .

On using (4.5) and only considering the passive component of the stress, relevant for determining the viscosity at large $Pe$ , one obtains

(4.6) $$\begin{eqnarray}\displaystyle \langle \unicode[STIX]{x1D70E}_{E}^{b}\rangle _{33}^{p}|_{outer} & = & \displaystyle \frac{\unicode[STIX]{x03C0}(nL^{3})}{18\ln \unicode[STIX]{x1D705}}\left[\frac{(\exp (3t)+2)(8\exp (3t)+1)}{4(\exp (3t)-1)^{2}}\right.\nonumber\\ \displaystyle & & \displaystyle -\left.3\exp (3t)\frac{(5\exp (3t)+4)\tan ^{-1}[(\exp (3t)-1)^{1/2}]}{4(\exp (3t)-1)^{5/2}}\right].\end{eqnarray}$$

In figure 8(a), we plot the extensional viscosity obtained numerically in § 2.2 along with (4.6), for $Pe_{r}=100$ , 200 and 300, congruent with the discussion in § 3.1. The analytical and numerical evolutions are largely in agreement. The numerical plateaus lie slightly below the analytical one owing to the effects of rotary diffusion, which eventually arrests the convective evolution, yielding an $\unicode[STIX]{x1D6FA}$ that is a narrow Gaussian with an $O(Pe_{r}^{-1/3})$ angular spread about the extensional axis, instead of the delta function.

Figure 8.

Comparison of the numerical and analytical high- $Pe$ relaxations. (a) Axisymmetric extensional flow, (b) simple shear flow at $Pe_{r}=100$ , 200 and 300. The solid line represent the analytically obtained stress, equations (4.6) and (4.10), respectively, while the dashed, dotted and dash-dotted lines represent the numerical solutions.

4.1.2 Simple shear flow

We choose a slightly different coordinate system from that in § 2.1.2. In the analysis below, the polar angle $\unicode[STIX]{x1D703}$ is measured from the flow axis and the azimuthal angle $\unicode[STIX]{x1D719}$ from the gradient axis in the gradient–vorticity plane. The choice is motivated by the meridional nature of the infinite-aspect-ratio Jeffery trajectories. In the chosen coordinates, $\unicode[STIX]{x1D719}$ only acts to parametrize the trajectories, with only $\unicode[STIX]{x1D703}$ varying along a given trajectory. Using $\dot{\boldsymbol{p}}=-\sin ^{2}\unicode[STIX]{x1D703}\cos \unicode[STIX]{x1D719}\,\hat{\unicode[STIX]{x1D73D}}$ , equation (4.2), now takes the form:

(4.7) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\,\unicode[STIX]{x1D6FA}_{o}}{\unicode[STIX]{x2202}t}-\sin ^{2}\unicode[STIX]{x1D703}\cos \unicode[STIX]{x1D719}\frac{\unicode[STIX]{x2202}\,\unicode[STIX]{x1D6FA}_{o}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}-3\sin \unicode[STIX]{x1D703}\cos \unicode[STIX]{x1D703}\cos \unicode[STIX]{x1D719}\,\unicode[STIX]{x1D6FA}_{o}=0. & & \displaystyle\end{eqnarray}$$

Again, use of the method of characteristics gives:

(4.8a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\unicode[STIX]{x1D703}}{\text{d}t}=-\sin ^{2}\unicode[STIX]{x1D703}\cos \unicode[STIX]{x1D719}, & \displaystyle\end{eqnarray}$$

(4.8b ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\unicode[STIX]{x1D6FA}_{o}}{\text{d}t}=3\sin \unicode[STIX]{x1D703}\cos \unicode[STIX]{x1D703}\cos \unicode[STIX]{x1D719}\unicode[STIX]{x1D6FA}_{o}. & \displaystyle\end{eqnarray}$$

The solution to (4.8a ), gives a characteristic curve of the form $\cot \unicode[STIX]{x1D703}=t\cos \unicode[STIX]{x1D719}+\unicode[STIX]{x1D703}^{\prime }$ , where $\unicode[STIX]{x1D703}^{\prime }$ is a constant of the integration. Using this to solve (4.8b ), we get:

(4.9) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FA}_{o}=\frac{1}{4\unicode[STIX]{x03C0}}\left(\frac{1}{1+t^{2}\cos ^{2}\unicode[STIX]{x1D719}\sin ^{2}\unicode[STIX]{x1D703}-2t\cos \unicode[STIX]{x1D719}\sin \unicode[STIX]{x1D703}\cos \unicode[STIX]{x1D703}}\right)^{3/2}. & & \displaystyle\end{eqnarray}$$

Finally, the shear viscosity determined using (4.9), is

(4.10) $$\begin{eqnarray}\displaystyle \langle \unicode[STIX]{x1D70E}_{SS}^{b}\rangle _{12}^{p}=\frac{\unicode[STIX]{x03C0}(nL^{3})}{6\ln \unicode[STIX]{x1D705}}\int \text{d}\boldsymbol{p}\sin ^{2}\unicode[STIX]{x1D703}\cos ^{2}\unicode[STIX]{x1D703}\cos ^{2}\unicode[STIX]{x1D719}\unicode[STIX]{x1D6FA}_{o}, & & \displaystyle\end{eqnarray}$$

the integral above being evaluated numerically. In figure 8(b), we plot the passive viscosity obtained numerically in § 2.2, for $Pe_{r}=100$ , 200 and 300, along with (4.10). Here, the initial peak in the viscosity is convective in origin, and as explained in § 3, is associated with the transient accumulation of swimmer orientations along the extensional axis on time scales of $O(\dot{\unicode[STIX]{x1D6FE}}^{-1})$ . The viscosity predicted by (4.10), captures this peak, and therefore, is in good agreement with the large- $Pe$ results. However, in contrast to the numerical solution which exhibits an oscillatory approach towards a finite $Pe$ -dependent plateau, equation (4.10) monotonically decays to zero, because a perfectly flow-aligned orientation distribution at steady state leads to a vanishing viscosity. The oscillations are a signature of rotary diffusion becoming important at large times, as is evident from the lack of collapse of the stress oscillations for different $Pe$ on a convective time scale.