2 Stochastic hydrodynamic model

Consider a passive particle submerged in an inviscid fluid in two dimensions. A classical problem in hydrodynamics is the potential flow past a cylinder moving along a straight line and the explicit formula for the 2-D streamfunction is available (Maxwell Reference Maxwell1869). Consequently, the drift experienced by the particle can be readily computed by integrating the velocities in time. It was shown by Thiffeault & Childress (Reference Thiffeault and Childress2010) that the total particle displacement due to infrequent encounters with a dilute suspension of cylinders swimming in random directions can be modelled by the linear superposition

(2.1) $$\begin{eqnarray}\boldsymbol{x}(t)=\boldsymbol{x}_{0}+\mathop{\sum }_{k=1}^{M(t)}{\it\Delta}_{{\it\lambda}}(a_{k},b_{k})\hat{\boldsymbol{r}}_{k},\end{eqnarray}$$

where $\boldsymbol{x}(t)$ is the particle displacement vector at time $t$ , $\boldsymbol{x}_{0}$ is its initial position, $(a_{k},b_{k})$ are the impact parameters imposed by the $k$ th encountered swimmer who moves for a fixed distance ${\it\lambda}$ in the random direction $\hat{\boldsymbol{r}}_{k}$ and $M(t)$ is the number of encounters as a function of time. With proper averaging, the effective diffusivity of the scalar field is

(2.2) $$\begin{eqnarray}{\it\kappa}:=\frac{\langle |\boldsymbol{x}(t)-\boldsymbol{x}_{0}|^{2}\rangle }{4t}=\frac{2Un}{{\it\lambda}}\int _{0}^{\infty }\int _{-\infty }^{\infty }{\it\Delta}_{{\it\lambda}}^{2}(a,b)\,\text{d}b\,\text{d}a,\end{eqnarray}$$

where $U$ is the constant speed of the cylinders and $n$ is the number density of the swimmers. To compute the individual drift ${\it\Delta}_{{\it\lambda}}$ , one only needs to differentiate the streamfunction of the potential flow past a cylinder and integrate the velocities in time. When the swimming distance ${\it\lambda}$ is much larger than the cylinder size $\ell$ and can be assumed to be infinite, Thiffeault & Childress (Reference Thiffeault and Childress2010) obtained

(2.3) $$\begin{eqnarray}{\it\kappa}_{s}\approx 1.19Un\ell ^{3}.\end{eqnarray}$$

The accuracy of the approximation (2.1) relies on the dilute assumption: the number density $n$ has to be small so that the interaction between the cylinders in the potential flow is negligible. In other words, the drift imposed upon the passive particle at any instance of time comes predominantly from one swimmer. This is violated when schools of multiple cylinders that are close to each other exist. In this paper, we look at the simplest non-trivial scenario of schooling: the swimmers form dilute, well-separated schools while within each school, two identical cylinders stay close to each other and move synchronously with identical speed, duration and direction. A diagram of each encounter between the passive particle and a school pair is illustrated in figure 1. This is very similar to Thiffeault & Childress (Reference Thiffeault and Childress2010) and Lin et al. (Reference Lin, Thiffeault and Childress2011) with two more parameters introduced: the separation between two cylinders, $2L$ ( $L\geqslant \ell$ ), and the inclination angle between the swimming direction and the line connecting the cylinder centres, ${\it\theta}$ .

Figure 1. The drift displacement ${\it\Delta}={\it\Delta}_{{\it\lambda}}(a,b,L,{\it\theta})$ induced by a pair of schooling cylinders. The impact parameters $a$ and $b$ denote the initial perpendicular distance and the horizontal distance, respectively, from the start of the particle trajectory (filled dot) to the midpoint between the cylinders. The configuration parameter $L$ is defined as half of the separation between the cylinder centres whereas ${\it\theta}$ is the inclination angle from the swimming direction to the line connecting the centres. The particle stops at the hollow dot when two cylinders finish their gliding of length ${\it\lambda}$ with constant speed $U$ .

3 Method of image doublets for potential flow past two cylinders

To extend the work of Thiffeault and Childress we now derive the streamfunction for the potential flow generated by two cylinders moving synchronously, as shown in figure 1. It is easy to see that a simple superposition of two doublets would distort the impermeable boundaries from being circular, which inspired several methods to compensate for the inter-cylinder effects, including conformal mapping (Crowdy Reference Crowdy2006), elliptic function theory (Johnson & McDonald Reference Johnson and McDonald2004) and the method of image doublets (Carpenter Reference Carpenter1958; Dalton & Helfinstine Reference Dalton and Helfinstine1971). Here we demonstrate the method of image doublets due to its simplicity and derive the streamfunction and velocities from complex analysis.

The basic idea is to construct an infinite series of image doublets with decreasing strength for each cylinder. Within each series, the first, zeroth-order doublet represents the unperturbed cylinder and the $k$ th ( $k\geqslant 1$ ) image offsets the boundary distortion caused by $(k-1)$ th-order image in the other series. Finally, the total complex potential is simply the sum of a uniform flow and two series

(3.1) $$\begin{eqnarray}w={\it\phi}+\text{i}{\it\psi}=-Uz+\mathop{\sum }_{k=0}^{\infty }w_{1,k}+\mathop{\sum }_{k=0}^{\infty }w_{2,k}.\end{eqnarray}$$

Here the uniform flow at infinity is moving from right to left and the convergence of these series is guaranteed by the decay of the doublet strength in each series.

Next we derive the formulas for the doublets $w_{1,k}$ and $w_{2,k},k=0,1,\ldots .$ This is equivalent to the determination of the position and the strength for each doublet. Without loss of generality, for each encounter illustrated in figure 1 we set up a co-moving, complex $z$ -plane with the midpoint between the two cylinders being the origin, and with the swimming direction being the positive real ( $x$ ) axis.

It is known that for the zeroth-order doublets that model the potential flow past a cylinder (Acheson Reference Acheson1990)

(3.2) $$\begin{eqnarray}w_{j,0}=-\frac{U\ell ^{2}}{z-z_{j,0}},\quad j=1,2,\end{eqnarray}$$

with $z_{j,0}=(-1)^{\,j-1}L\,\text{e}^{\text{i}{\it\theta}}=\pm L(\cos {\it\theta}+\text{i}\sin {\it\theta})$ the coordinates of the cylinder centres that are symmetric with respect to the origin moving with the cylinders. For the next order, since the image doublet $w_{1,1}$ offsets the circular boundary perturbation induced by doublet $w_{2,0}$ (within cylinder 2) around doublet $w_{1,0}$ (within cylinder 1), the following restriction should be imposed on the imaginary parts of the doublet potentials:

(3.3) $$\begin{eqnarray}\text{Im}(w_{1,1}+w_{2,0})=-U\ell ^{2}\text{Im}\left(\frac{s_{1,1}}{z-z_{1,1}}+\frac{1}{z-z_{2,0}}\right)=\text{constant}\end{eqnarray}$$

in which $s_{1,1}$ and $z_{1,1}$ , the relative strength and position of the image doublet respectively, are chosen as follows (Carpenter Reference Carpenter1958):

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle s_{1,1}=-\frac{\ell ^{2}\text{e}^{2\text{i}{\it\theta}}}{|z_{1,0}-z_{2,0}|^{2}}=-\frac{\ell ^{2}\text{e}^{2\text{i}{\it\theta}}}{4L^{2}}, & \displaystyle\end{eqnarray}$$

(3.5) $$\begin{eqnarray}\displaystyle & \displaystyle z_{1,1}=z_{1,0}+\frac{\ell ^{2}}{\overline{z_{2,0}-z_{1,0}}}=\text{e}^{\text{i}{\it\theta}}\left(L-\frac{\ell ^{2}}{2L}\right), & \displaystyle\end{eqnarray}$$

where $\overline{\;\cdot \;}$ denotes complex conjugacy. Note how the strength of the image decays according to an inverse square law for variable $L$ and how it lies on the line connecting two cylinder centres. Similarly, the formula for $w_{2,1}$ , the image doublet that restores the boundary distortion around cylinder 2 by $w_{1,0}$ can be derived.

To summarise, the first-order image doublets in the complex potential (3.1) are

(3.6) $$\begin{eqnarray}w_{j,1}=-\frac{U\ell ^{4}\text{e}^{2\text{i}{\it\theta}}}{4L^{2}\left[z+(-1)^{j}\text{e}^{\text{i}{\it\theta}}\left(L-{\displaystyle \frac{\ell ^{2}}{2L}}\right)\right]},\quad j=1,2.\end{eqnarray}$$

In fact, all higher-order image doublets can be derived as above to balance corresponding lower-order images in the same inductive fashion but with more tedious details. However, as we will see in the next sections, first- and second-order images are sufficient for the purpose of computing effective diffusivity. Furthermore, this procedure can be readily generalised for cylinders of different sizes and for more than two cylinders (Dalton & Helfinstine Reference Dalton and Helfinstine1971). We elect to postpone the discussion for these scenarios to future work since preliminary results show that more complicated configurations do not lead to significantly new phenomena in current context.

It should also be noted that the requirement (3.3) only enforces that the circular boundaries of the cylinders are impermeable streamlines. The constant on the right-hand side is generally non-zero and therefore the cylinders may be subject to lift and drag forces. For more detailed derivations readers can refer to Dalton & Helfinstine (Reference Dalton and Helfinstine1971). To study the schooling effects, here we assume that there are ‘internal’ forces constantly exerted on the cylinder pair to maintain their relatively stationary positions.

4 Results for the effective diffusivity

With the complex potential formulas presented above, we now compute the effective diffusivity ${\it\kappa}$ defined in (2.2) for schools of two cylinders with two configuration parameters: cylinder separation $L$ and inclination ${\it\theta}$ . Figure 2 summarises the results for three typical formations: ${\it\theta}=0$ (‘chasing’), ${\rm\pi}/4$ (‘tilting’) and ${\rm\pi}/2$ (‘sweeping’). The horizontal axis is the distance between cylinder centres normalised by the cylinder radius; the vertical axis is the effective diffusivity normalised by twice the reference value (2.3). The factor $2$ is introduced to highlight the nonlinearity in the schooling enhancement compared with simply doubling the swimmer density in the original model. These curves are numerically generated by truncating each of the two series in the potential (3.1) to three terms and with constants $U=\ell =1$ , $n=10^{-3}$ and ${\it\lambda}=100$ . It has been verified for a wide range of parameter settings that including more terms can change the result no more than $1\,\%$ due to the fast decay of image doublets. In fact, keeping only two terms in each series, namely, considering only the zeroth-order doublets $w_{j,0}$ and their first-order corrections $w_{j,1}$ , $j=1,2$ , recovers more than $96\,\%$ of the effective diffusivity. Although the truncation does result in slight distortions in the cylindrical boundaries of the swimmers, and thus some error in computing the particle drift very close to the boundaries, their contribution to the integrated effective diffusivity is negligible.

Figure 2. Normalised effective diffusivity ${\it\kappa}$ as a function of cylinder separation $L$ . For different in-school formations, the dependence of the effective diffusivity on the separation parameter exhibits opposite behaviours: when ${\it\theta}=0$ , the diffusivity increases as the distance between two cylinders grows while for ${\it\theta}={\rm\pi}/4$ and ${\it\theta}={\rm\pi}/2$ the monotonicity is reversed. Here the number density $n=10^{-3}$ , swimming speed $U=1$ and the cylinder radius $\ell =1$ are all kept constant.

Here we observe two opposite behaviours of the effective diffusivity as a function of the separation parameter: for the chasing case ( ${\it\theta}=0$ ), ${\it\kappa}$ is a strictly increasing function of $L$ while it is a strictly decreasing function in tilting ( ${\it\theta}={\rm\pi}/4$ ) and sweeping ( ${\it\theta}={\rm\pi}/2$ ) formations. Moreover, sweeping schools yield a much bigger boost than tilting ones when $L$ is small. In all three cases, the dependence is nonlinear in that ${\it\kappa}$ varies rapidly for $L/\ell <2$ and approaches an asymptotic, constant value as $L$ gets large.

A straightforward intuition can be applied to explain the two limiting cases for large $L$ : When a school of two swimmers chase through the fluid with enough separation, each encounter with the particle can be well approximated by two sequential kicks by the leading then the trailing cylinder from exactly the same direction. This long-range superposition doubles combined particle drift which implies a 4-fold increase in the squared displacement ${\it\Delta}^{2}$ and consequently ${\it\kappa}\approx 4{\it\kappa}_{s}$ . On the other hand, in the tilting and sweeping cases where ${\it\theta}$ is significantly different from 0, the two cylinders no longer cooperate with each other when they are far apart and they act on the particle independently. The resulting effective diffusivity is then simply a linear extrapolation of the value from independent swimmers and therefore ${\it\kappa}\approx 2{\it\kappa}_{s}$ . In other words, the swimmers achieve no extra advantage by schooling together under these configurations. Of course, this asymptotic analysis is only valid when $L$ is still small compared to the average distance between two schools so the dilute assumption is still accurate.

More careful investigation is required when $\ell <L<2\ell$ . For the chasing configuration, figure 3 illustrates how the particle drift depends on the separation $L$ . When the two cylinders are slightly apart ( $L/\ell =1.1$ ), the overlaid particle trajectory shows that the combined drift is only approximately equal to the single-cylinder value. In this case, the schooling in fact suppresses the mixing efficiency since doubling the number of swimmers by schooling does not double the effective diffusivity. This regime can also be identified in figure 2 in which part of the solid curve ( ${\it\theta}=0$ ) falls below unit value, namely, ${\it\kappa}/2{\it\kappa}_{s}<1$ . As $L$ grows, the particle trajectory essentially becomes a superposition of two successive encounters as mentioned before and shown by the other overlaid trajectory in the figure ( $L/\ell =2$ ). Therefore, schools in chasing formations can only enhance mixing significantly when the in-school separation is large enough.

Figure 3. The effect of cylinder separation $L$ on the drift displacement ${\it\Delta}$ when ${\it\theta}=0$ . Here the impact parameters $a/\ell =0.1$ and $b/\ell =5$ are fixed. Two typical trajectories are overlaid. When $L$ is small, two cylinders operates nearly as one for the particle; when $L$ is large enough, two cylinders induce a drift displacement that almost doubles the single-cylinder value.

5 Amplified active region and short-range coupling for ${\it\theta}\neq 0$

Finally, we investigate the source of the nonlinearly enhanced diffusivity when ${\it\theta}\neq 0$ . In this section we fix $L=1.2$ without loss of generality. First we examine the dominating contribution of the transformed ${\it\kappa}$ integral in the $\log (a/\ell ){-}(b/{\it\lambda})$ parameter plane

(5.1) $$\begin{eqnarray}{\it\lambda}^{-1}\ell ^{-3}\int _{-\infty }^{\infty }\int _{0}^{\infty }{\it\Delta}_{{\it\lambda}}^{2}(a,b)\,\text{d}a\,\text{d}b=\iint _{\mathbb{ R}^{2}}\ell ^{-3}a{\it\Delta}_{{\it\lambda}}^{2}(a,b)\,\text{d}\log (a/\ell )\,\text{d}(b/{\it\lambda})\end{eqnarray}$$

by comparing in figure 4 the distributions of the dimensionless integrand $\ell ^{-3}a{\it\Delta}_{{\it\lambda}}^{2}(a,b)$ (Lin et al. Reference Lin, Thiffeault and Childress2011) for the non-schooling case and for the three schooling formations. The active region in each panel is the area where the integrand value is significantly non-zero (shown in darker colour), or equivalently, is the set of $(a,b)$ such that $\ell ^{-3}a{\it\Delta}_{{\it\lambda}}^{2}(a,b)\geqslant {\it\gamma}$ where the significance ${\it\gamma}$ is set to be $0.3$ here. It should be noted that any reasonable choice of ${\it\gamma}$ supports the subsequent arguments.

Figure 4. Distributions of integrand $\ell ^{-3}a{\it\Delta}_{{\it\lambda}}^{2}(a,b)$ in (5.1) for (a) a single cylinder, and for a two-cylinder school in three different configurations: (b) ${\it\theta}=0$ , (c)  ${\it\theta}={\rm\pi}/4$ and (d) ${\it\theta}={\rm\pi}/2$ . The significant values (plotted with darker shading) in cases (a) and (b) spread over $a\ll \ell ,b/{\it\lambda}\in (0,1)$ while they are greatly intensified and localised near $a\approx \ell$ in cases (c) and (d). In all figures $U=\ell =1,{\it\lambda}=100$ and for all schooling configurations $L/\ell =1.2$ .

We can see that the single-cylinder case and the chasing case (figure 4 a,b) share a strong similarity. In both cases the integral is dominated by the active region $\log (a/\ell )\lesssim -0.6,b/{\it\lambda}\in (0,1)$ which corresponds to the ‘head-on’ collisions near $a=L\sin {\it\theta}=0$ (see figure 1 with ${\it\theta}=0$ ) in physical coordinates. And a chasing school produces longer drifts and thus a larger diffusivity. By contrast, for tilting and sweeping schools (figure 4 c,d) the dominating contributions come from a shifted and more localised region near $\log (a/\ell )=0$ , or $a=L\sin {\it\theta}=O(1)$ . As it moves away from $a\ll \ell$ , the intensity in this region is greatly amplified to an extent that the integrated diffusivity achieves a nonlinear growth for ${\it\theta}\neq 0,\ell <L<2\ell$ as shown in figure 2. In fact, the integrand reaches its maximal value of 10 (or 16) when ${\it\theta}={\rm\pi}/4$ (or ${\it\theta}={\rm\pi}/2$ ) for typical values of $b$ while the maximum is only 0.6 in the non-interacting case. This magnification factor is much more than enough to compensate for the active area localisation when integrating ${\it\kappa}$ in (5.1) although, interestingly, this region actually expands in the dimensional $a$ – $b$ plane which will be mentioned again. Both the magnification in strength and the expansion in dimensional space are much more significant for ${\it\theta}={\rm\pi}/2$ than for ${\it\theta}={\rm\pi}/4$ and thus the advantage of sweeping over tilting entails.

The physical origin of the above analysis is attributed to the short-range coupling via which tilting and sweeping schools with small $L$ produce large displacements and super-linearly enhanced diffusivity. This coupling involves the simultaneously impacts on the particle from both cylinders when ${\it\theta}\neq 0$ which are in the same order of magnitude. Furthermore, the coupling effect intensifies as ${\it\theta}$ increases towards ${\rm\pi}/2$ , or as $L$ decreases which agrees with what we have seen in figure 2. However, this effect is generally negligible when ${\it\theta}=0$ since the subdominant, ‘farther’ cylinder only provides a contribution of the order of $O(1/(1+2L/\ell )^{2})$ of what the dominant, ‘closer’ one enforces.

Figure 5 visualises this mechanism by superimposing typical trajectories for different values of ${\it\theta}$ and $a$ . The key observation one can make here is that the two cylinders within each of the ${\it\theta}\neq 0$ schools (centre and right panels) reinforce each other not just by invoking long drifts but also by slowing down the drift decay away from the head-on positions, $a=L\sin {\it\theta}$ , such that particles with a broader range of impact parameters experience long drifts in comparison with the chasing case as noted in the previous paragraph. It should also be noted that these enhanced drifts deviate from the close loops formed with same impact parameters under non-interacting settings (Lin et al. Reference Lin, Thiffeault and Childress2011; Pushkin et al. Reference Pushkin, Shum and Yeomans2013).

Figure 5. Typical trajectories of a target particle in the fixed laboratory frame for different values of $a$ due to a two-cylinder school moving from left to right with $U=\ell =1$ , ${\it\lambda}=100$ and $b/{\it\lambda}=0.3$ . The initial positions of the particle are marked by solid dots and the final positions by hollow dots. In each panel, solid dots from top to bottom correspond to $a/\ell =1.2,0.9,0.6,0.3,0.05$ . The cylinder separation is fixed in all three configurations at $L/\ell =1.2$ . Note how the fore–aft symmetry of a trajectory is completely broken by a tilting school as seen in the centre panel: (a) ${\it\theta}=0$ , (b) ${\rm\pi}/4$ , (c) ${\rm\pi}/2$ .

6 Schooling effects of three cylinders

Our next step is to explore the mixing effects of schools of three or more cylinders. This is particularly natural in light of our motivating problem, ocean biogenic mixing, since in reality a fish school contains up to thousands of individuals. In theory, the method of image doublets can be applied to arbitrary number of cylinders in two dimensions. For example, for a school of three cylinders the construction of image doublet series analogous to § 3 would start from three zeroth-order doublets, to six first-order image doublets (two first-order images within each zeroth-order doublet to balance the other two zeroth-order ones), to twelve second-order images and so on.

We have therefore conducted a preliminary study on triple-cylinder schools as extensions to two-cylinder formations. It is to be expected that the effective diffusivity would have more subtle dependence on the parameters due to the extra degrees of freedom characterising in-school configuration. Again, we truncated the infinite series to first-order images and more than 95 % of the mixing effects was captured in all of our test cases. In particular, in figure 6 we show three representative configurations for an equilateral, three-cylinder school with side lengths $2L=2.4\ell$ and these configurations all produce super-linear mixing enhancement ( ${\it\kappa}/3{\it\kappa}_{s}>1$ ). In fact, the large- $L$ limit of the normalised effective diffusivity, ${\it\kappa}/3{\it\kappa}_{s}$ , can go up to 3 in contrast to 2 in the two-cylinder cases shown in figure 2. Not surprisingly, this can be achieved by a chasing formation in which three cylinders align and move along a straight line just as the two-cylinder case. This is demonstrated by the solid curve in figure 7 which plots the effective diffusivity as a function of separation parameter $L$ under two types of formation shown as inlets. The other curve (dashed) is produced by a triangular formation with an evident sweeping feature. In summary, these results show striking similarity with their counterparts in figure 2 and thus add greatly to our confidence that the two mechanisms identified in the study for two-cylinder cases remain crucial in more crowded schools.

Figure 6. Three configurations for a three-cylinder school. In each configuration the cylinders form an equilateral triangle with $L$ fixed at $1.2\ell$ and three cases differ in the triangle orientation relative to the moving direction: (a) ${\it\kappa}/3_{{\it\kappa}_{s}}=2.36$ , (b) 2.36, (c) 2.52.

Regardless of whether it is plausible to extrapolate this into a quadratic rule, namely,

(6.1) $$\begin{eqnarray}\max \frac{{\it\kappa}}{{\it\kappa}_{s}}=N^{2}\end{eqnarray}$$

over all possible in-school configurations where $N$ is the school size, or whether the chasing formation is always the optimal, a more detailed study would be required. As suggestive as figures 2 and 7 may seem, we are yet to conclude that the maximal diffusivity is indeed 3 and that it is uniquely realised by a chasing formation attributing to the difficulty in exhausting and parametrising all possible three-cylinder configurations, such as triangles with unequal side lengths, not to mention the resulting formulas for image doublets that are much more tedious and less tractable. Complications of this nature are further multiplied when we consider four or more cylinders within one school.

However, the above results and analyses do provide some insights, at least under the settings of a planar, potential flow, to how in-school interactions between cylinders can translate to enhanced mixing. In both two-cylinder and three-cylinder settings, the long-range superposition (successive kicks) and short-range coupling (amplified active region) have been shown in figure 7 to be effective mechanisms in terms of producing super-linear enhancements. With the inclusion of even more cylinders, we conjecture that these two would remain the dominant contributors although their simultaneous operation would definitely add to the nonlinearity and thus the complexity of the problem. For example, we compute the normalised diffusivity for the four-cylinder configuration illustrated in figure 8 and found that ${\it\kappa}/4{\it\kappa}_{s}\rightarrow 2.45$ as $L\rightarrow \infty$ which significantly deviates from unity seen for non-chasing two- or three-cylinder formations. Due to the scope of this manuscript, we elect to postpone further study along this direction to future work.

Figure 7. Normalised effective diffusivity as a function of cylinder separation $L$ for three-cylinder schools in two types of formation in which the behaviour of the effective diffusivity follows similar patterns as in figure 2.

Figure 8. Square formation for a four-cylinder school. The normalised diffusivity is found to approach 2.45 as $L$ grows.

7 Discussion and future work

In this work, we focus on the mixing efficiency of a schooling pair of two cylinders moving synchronously in a potential flow and conclude that chasing and sweeping schools are generally more advantageous than tilting ones. Furthermore, we illustrate the two physical mechanisms, long-range superposition ( ${\it\theta}=0$ ) and short-range coupling ( ${\it\theta}\neq 0$ ), through which these schools generate large drift displacements and thus a super-linear growth in the effective mixing diffusivity compared with non-interacting cylinders. Among others, we will improve our understanding to the impact of schooling in the context of ocean mixing and of other realistic scenarios in several aspects.

First is to study schooling effects in other fluid regimes and geometries, such as Stokes flow, a classical model for slow moving bodies in viscous fluids and other three-dimensional flows. We limited our discussion to two-dimensional schools due to the simplicity in the series representation of image doublets and in the physical explanation to the enhanced diffusion. Under axisymmetric scenarios, Lin et al. (Reference Lin, Thiffeault and Childress2011) considered the mixing effects of a dilute suspension of independent Stokesian squirmers (Lighthill Reference Lighthill1952; Blake Reference Blake1971) in three dimensions and established connections to various biofluids (Ishikawa & Pedley Reference Ishikawa and Pedley2007; Drescher et al. Reference Drescher, Leptos, Tuval, Ishikawa, Pedley and Goldstein2009; Leptos et al. Reference Leptos, Guasto, Gollub, Pesci and Goldstein2009; Guasto, Johnson & Gollub Reference Guasto, Johnson and Gollub2010), along with analogous results for independent spheres in a three-dimensional potential flow. More complicated boundary conditions and the loss of axisymmetry generally pose substantial challenges to extend the current theoretical framework. One way to explore schooling effects under those settings is to first study special, axisymmetric formations of three-dimensional schools for which one can still characterise the flow field with streamfunctions and singularities, although the construction of image doublets would be much more involved. As a qualitative conjecture, we expect that the mixing enhancement would be weaker in three dimensions, analogous to the comparison between cases in two and three dimensions documented in Thiffeault & Childress (Reference Thiffeault and Childress2010) due to the extra dimension and faster velocity decay. On the other hand, with viscosity introduced as a realistic fluid condition, we expect even higher mixing enhancement as a result of stronger coupling between swimmers.

Alternatively, one should also include more swimmers within each school in the context of studying ocean biogenic mixing, since in reality a fish school contains up to thousands of individuals. We have seen in § 6 how straightforward generalisations of the method of image doublets for three-cylinder schools yielded complicated formulas and incomplete results. For even larger schools, we may need to resort to a simple, macroscopic effective model to somehow parametrise the microscopic in-school interactions. Specifically, a school of small fishes can behave effectively like a large fish at ocean mixing which, unlike a small fish, exhibits considerable mixing efficiency (Visser Reference Visser2007). What we can reasonably conjecture for a synchronised school of thousands of individual swimmers is that under certain formations (relatively small swimmer separations perpendicular to the swimming direction as in sweeping two-cylinder schools, and not too small separation in the swimming direction as in chasing two cylinder schools), the two mechanisms enhancing mixing discussed before would coexist and cooperate in a way that demands further investigation. In fact, figure 7 has shown that three-cylinder schools that possess dominant chasing or sweeping components do indeed produce greater enhancements.

Additionally, it would be interesting to be able to identify the critical inclination ${\it\theta}^{\ast }>0$ , if it exists, at which the curves shown in figure 2 bifurcate from the increasing behaviour for ${\it\theta}=0$ to the decreasing behaviour for ${\it\theta}={\rm\pi}/4$ as the separation $L$ grows.