Correlation properties of collective motion in bacterial suspensions

The study of collective motion in bacterial suspensions has been of significant recent interest. To better understand the non-trivial spatio-temporal correlations emerging in the course of collective swimming in suspensions of motile bacteria, a simple model is employed: a bacterium is represented as a force dipole with size, through the use of a short-range repelling potential, and shape. The model emphasizes two fundamental mechanisms: dipolar hydrodynamic interactions and short-range bacterial collisions. Using direct particle simulations validated by a dedicated experiment, we show that changing the swimming speed or concentration alters the time scale of sustained collective motion, consistent with experiment. Also, the correlation length in the collective state is almost constant as concentration and swimming speed change even though increasing each greatly increases the input of energy to the system. We demonstrate that the particle shape is critical for the onset of collective effects. In addition, new experimental results are presented illustrating the onset of collective motion with an ultrasound technique. This work exemplifies the delicate balance between various physical mechanisms governing collective motion in bacterial suspensions and provides important insights into its mesoscopic nature.


Introduction
The phenomenon of collective motion has been an active area in the last few years. Experiments have demonstrated that in the collective state a bacterial suspension can exhibit remarkable properties such as enhanced diffusivity, the formation of sustained whorls and jets, and the ability to extract useful work [1][2][3][4][5]. Specifically in suspensions of rear actuated swimmers (socalled pushers) exemplified by bacteria, correlated mesoscopic collective motion is observed at length scales up to ten times larger than that of a given particle while particle velocities exceed that of an isolated swimmer [2,6].
While specific biological mechanisms such as chemotaxis and chemical signaling can be ruled out in most of the experiments on collective swimming in concentrated bacterial suspensions, the role of the physical mechanisms governing collective motion is still under discussion. A variety of theoretical models have been considered recently, with fundamentally different mechanisms attributed to the onset of collective motion, from flagella reversal dynamics [4], pure hydrodynamic [7][8][9][10][11] to pure collisional [4,[12][13][14][15][16] interactions. While the existing models were able to reproduce some key experimental observations, such as the onset of collective motion above a certain critical concentration [6,11,17], but the agreement with experiment remains incomplete and mostly qualitative.
A variety of computational approaches were suggested to model individual swimmers, such as slender bodies [9], swimming dumbbells [14][15][16], hard self-propelled rods [12], self-propelled hydrodynamic point dipoles [18], squirmers [19] or more computationally intensive 'composite' swimmers [20,21]. Each model, emphasizing different aspects of the balance between hydrodynamic forces and steric collisions, often yields to somewhat different conclusions on the critical concentration for the onset of collective motion, swimming speed in the collective state or its correlation properties. In this paper, by introducing a minimal yet realistic model of a swimmer, we investigate numerically the combined effect of hydrodynamic interactions and collisions on the collective dynamics of a bacterial suspension. A dedicated experimental study was conducted to determine the time required to establish the collective state from specific initial conditions (well-mixed or large vortex), and spectral and spatiotemporal velocity correlations in the corresponding collective state. The main goal of this study is to provide further understanding of the phenomena discovered in our recent experiments: the velocity correlation length and properly rescaled correlation time of the collective motion are essentially independent of the velocity and concentration of bacteria [22].
The structure of the paper is the following. First, we briefly present an overview of experimental techniques used to study the onset of collective motion. Then, we introduce the model for the suspension where a bacterium is represented as a point force dipole with size and shape, where the force of self-propulsion induced by the beating of flagella is balanced by the viscous drag force on the bacteria. The force dipole model for a bacterium has been experimentally verified by Drescher et al [13]. A novel element of our model is an additional realigning short-range torque turning bacteria away after the collision. Without this feature the bacteria form unphysical dense closely packed aggregates. Using this model we investigate the effects of the system parameters (swimming speed/propulsion force, concentration and tumbling rate) on the collective dynamics of the suspension. We show that our model is justified by capturing earlier experimental results on the transition from individual to collective motion by studying the effect of concentration on the correlation length and mean particle velocity [2,6]. Subsequently, we present the results of varying the system parameters in the simulation and compare with similar experimental results. Moreover, our theoretical modeling allows for the separation of the contributions of tumbling and swimming speed as well as the effect of the particle shape on the correlation length, which at present has not been done experimentally.

Brief description of experiment
Our experimental setup, based on a free-standing thin film suspension of Bacillus subtilis, is similar to that described in [2,22]. Collective motion is quantified by measuring the bacterial velocity field, then the associated correlation function is defined and the correlation length and time extracted (see section A.3 of the appendix for the definition of the correlation function). Specifically, collective motion is characterized by an increased correlation length (quantifying the size) and mean particle velocity [2]. In addition, recent experiments demonstrate the effects of system parameters such as concentration, tumbling rate and swimming speed on the collective dynamics of the suspension [22,23]. These observations were made on both a three-dimensional thin film and quasi-two-dimensional (2D) film using a combination of both particle tracking velocimetry imaging technique for low concentrations and the optical flow method or particle image velocimetry (PIV) for high concentrations. The experiments were performed on bacterial suspensions where the swimming and tumbling rate were both controlled by varying the ratio of oxygen and nitrogen in the chamber containing the suspension. The correlation between these two modes of motion results from the simple observation that when deprived of oxygen bacteria swim more slowly and sensing unfavorable conditions they begin to tumble more frequently. New experimental results, not reported earlier in [22], are shown in figures 1-4 and 6.
The main experimental observation made in [22,23] was that the correlation time of the bacterial flow field in the collective state scaled inversely with the swimming speed or concentration while the correlation length remained constant for negligible amounts of noise. Thus suggesting that changing the bacterial swimming speed or concentration is analogous to changing the overall time scale. This follows from the fact that the magnitude of a bacterium's induced flow and the number of collisions between bacteria are proportional to their swimming speed and concentration. Each quantity is also proportional to the rate at which energy is injected into the system.
In a new experiment presented here we investigated the onset of collective swimming in a suspension of bacteria from an initially prepared isotropic disordered state. A small drop of bacterial suspension was placed on a microscope glass slide. An ultrasonic transducer was glued to the glass with epoxy, figure 1(a). A square wave frequency modulated signal (sweep) in the range 35-143 kH was applied to the transducer generating capillary waves on the open surface of the drop, with the typical wavelength of the order 60 µm, figure 1(b). The surface capillary waves effectively mix the suspension and completely destroy local collective motion (however, they also create a large scale rectified vortex, see the discussion below). We found that after the signal is turned off the bacterial swimming pattern is different from the stationary collective state and bacteria swim much slower (see movie 1 in the supplementary information (available from stacks.iop.org/NJP/15/105021/mmedia)). However, the swimming speed increases in time during the first 5 s, figure 4, and collective motion appears in ≈ 1 s, figure 1(d).
For detailed analysis of the onset of collective behavior we monitored the evolution of the power spectrum distribution P(k n ) = |F(k n )| 2 , where F(k n ) are the coefficients of the Fourier transformation of the velocity field and k n is the dimensionless wavenumber defined as k n =k n l, assuming l ≈ 10 −2 mm is an effective (bacterium body and flagella) bacterial length andk n is the dimensional wavenumber, both in experiment and simulations. Some results are presented in figure 2. We have found that while in the steady state, after the transition, the power spectra in the experiment and simulations are qualitatively similar, the details of the time evolution are different due to the specifics of experiment. We observed that the surface waves, excited by the US transducer, generate a rectified flow with a scale of a few millimeters. Thus, a large vortex with a big correlation length will be present initially, then slowly decreasing to the stationary value. These large scale flows decay in time very slowly due to the reduction of the effective viscosity of the bacterial suspension, [24]. To exclude the flow generated by the surface waves, we subtracted the average velocity in the microscope field of view from the velocity field. According to our experimental observation in figures 2(a) and (b), the average wavenumber is slightly increasing in time during the development of collective behavior. It can be evidence of the slow decay of a large scale flow generated by the surface waves to a flow with scales typical for collective swimming of the bacteria. In the simulations started from an initial vortex (vortex-like configuration of bacteria), figures 2(c) and (d), we see qualitatively the same behavior as in experiment; however, in simulations started from a random configuration of bacteria with no external large-scale flow we found that the average wavenumber is decreasing in time, figures 2(c) and (d).
In the steady state, both experiment and simulations yield the power spectra with a well-pronounced peak at a 'mesoscale' wavenumber k ≈ 0.3, compare figures 2(b) and (d).
We also investigated the energy spectra, E(k) ∼ k P(k), and obtained qualitatively similar results as in [12]. However, we do not think that it is appropriate to interpret such 0 min max 50 µm energy spectra as evidence for the intermediate power-like behavior and, correspondingly, 'mesoscale turbulence'-the interval in k-space for the power-like behavior is exceedingly small.
Figures 3(a) and (b) illustrate the velocity field and vorticity observed in experiment at two different swimming speeds for a concentration in the collective state to show the onset of vortex-like structures in higher volume fraction suspensions. One sees that for a small oxygen concentration (i.e. small swimming speed) the flow pattern is rather disordered due to random tumbling, with the characteristic size of the flow is on the order of the bacterial length. In contrast, for larger oxygen concentration (larger swimming speeds) the bacteria develop mesoscale patterns exhibiting collective motion manifested by self-organized vortices and jets with the characteristic size of about 6-10 bacterial lengths. Qualitatively similar behavior was also observed in our simulations, figures 3(c) and (d). Correspondingly, movies 1 and 2 in the supplementary information (available from stacks.iop.org/NJP/15/105021/mmedia) taken from simulation illustrate the visual similarity between the simulations and experiment for the dynamics of a suspension for both low and high swimming speeds.

Model
Here we introduce a simple computationally effective model to investigate the role of collisions and hydrodynamic interactions on the collective behavior of bacterial suspensions in two and three dimensions. Our model does not include buoyancy and oxygentaxis, and, therefore, cannot capture complex behavior such as bio-convection associated with collective motion in threedimensional (3D) suspensions [25,26]. Since bacteria can come in diverse shapes and sizes, using our simulations we will determine if the shape/aspect ratio effects change the collective motion. Additionally, bacteria do not swim in a straight line even in the absence of tumbling (the rotation of the flagella induces counter-rotation of the bacterial body), but this effect will be neglected for the purposes of the simulations.
In our previous theoretical work [18], we introduced a simple 3D model for semidilute bacterial suspensions which captured the qualitative effects of collective swimming on the effective viscosity observed in experiment. Simulations of this model performed on graphic processing units demonstrated, in agreement with experiment [24], an initial decrease and then subsequent increase in the effective viscosity with the increase in the concentration of bacteria over the full range of experimentally observed concentrations. In this paper we will use an extension of this model in two and three dimensions to study the dependence of the correlation length of collective motion as function of swimming speed, concentration, dipolar strength and geometrical confinement (two versus three dimensions).
We focus here on a semidilute suspension-the hydrodynamic interactions on each bacterium are represented as the sum of pairwise interactions. This approximation is exact for bacteria modeled as point dipoles. We consider two basic forms of interactions: bacterium with fluid (hydrodynamic) and inter-bacterial short-range interactions (soft collisions). The combination of hydrodynamic interactions and collisions has been shown to be crucial for collective motion [4,17]. The intuitive theoretical picture is that collisions align the particles and once aligned hydrodynamic interactions are enhanced. The main goal for the model is to capture properties of collective motion observed in suspensions of Bacillus subtilis, a pusher. This point dipole model can be extended to study suspensions of 'pullers' (swimmers propelled from the front) as well, such as the algae Chlamydomonas, by reversing the orientations of the point forces so that they point inward.
Due to the typical effective size (body and flagellum) of a bacterium (l = 10 µm) and isolated swimming speed (V 0 = 20 µm s −1 ), the Reynolds number Re = lV 0 /ν 1 is exceedingly small, where ν is the kinematic viscosity. Since a steady state flow is established on a timescale much smaller than the characteristic time scale of bacterial motion, the fluid is modeled as an incompressible steady Stokes flow neglecting inertial effects.
Two models will be proposed to handle the cases of a thin film quasi-2D suspension and a full 3D model. The first is designed to capture the thin film geometry representative of the experimental setup of Sokolov et al [2] and more recent experiment in thin microfluidic chamber by Dunkel et al [23]. The domain is a thin film of length and width L l and thickness w ∼ l approximately the characteristic length of a swimmer. While the experiments [2] were performed in a free-standing liquid film, the accumulation of bacteria metabolism products results effectively in no-slip boundary conditions for the fluid on the top and bottom surfaces of the film. In [23] the experiment is performed in suspension confined between two solid surfaces implying no-slip conditions as well. A suspension consists of N neutrally buoyant particles. Each particle will have two variables identifying it: the center of mass x = (x, y, z) and the orientation d = (d x , d y , d z ) directed along a given particle's axis of symmetry. We start from a confined thin film model for the fluid studied in [27].
We define the flow at the middle 2D plane of the thin film (z = 0) due to a point dipole bacterium with orientation d at the origin (|d| = 1) where r = |x| is the distance from the origin to the dipole, U 0 = ζ V 0 l 2 is the dipole strength where ζ is a constant determined by the shape/aspect ratio. This formula for a point dipole is derived from that of two oppositely oriented force monopoles in the limit as the distance between them goes to zero (cf [27]). Once obtained we find that the in-plane fluid velocity u x (x), u y (x) ∼ r −3 . Since the film has thickness approximately the size of a particle in the vertical dimension, we consider the suspension in two dimensions (x, y) for numerical simulations. We also consider a full 3D model for the suspension where the fluid is again modeled using the incompressible Stokes equation. The flow induced by a single bacterium is taken to be wherer = r/r , r = |r| is the distance vector relative to the center of dipole and U 0 = ζ V 0 l 2 . This scenario is representative of experiments done in the bulk of a suspension far from the boundary as in [6,12,22,23]. Particle motion is governed by the same equations of motion with (1) replaced by (2). The differences between the flows in the quasi-2D and full 3D are further explained in section A.2 of the appendix.
The fluid velocity for the suspension, u(x), depends on the N particle positions and orientations. The corresponding equations of motion, derived from the balance of forces and torques on each particle, are defined in section A.1 of the appendix. Due to the linearity of the Stokes equation, the flow due to N particles will be the sum of flows due to a single particle (1) or (2).
We study a full range of concentrations where we define the area fraction = π Nl 2 /4L 2 in 2D and the volume fraction = π Nl 3 /6L 3 in 3D which are the associated area/volume fractions for a sphere of radius l/2. Since we operate with point dipoles, the particle size, l, and the excluded volume are defined via the radius of the short-range generic repulsion potential σ LJ . In our simulations we used for convenience the truncated Lennard-Jones potential W ∼ [(σ LJ /r ) 12 − (σ LJ /r ) 6 ]. While the specific choice of the repulsive potential is not important (e.g. the exponential Yukawa potential instead of the Lennard-Jones potential was employed in [12]), the use of the truncated Lennard-Jones potential is computationally efficient 5 . We formally work with spherical particles for computational efficiency, but the shape of the particle is implicitly defined by the strength of the hydrodynamic dipole in (1)- (2). In addition, a certain realigning torque was added to ensure that the particles swim off after a collision, see sections A.2 and figure A.1 of the appendix. This effect is also related to the shape of particle.
Simulations were run for N identical particles, N ∼ 10 4 -10 5 , where a typical simulation contains between 15 000 and 40 000 particles depending on the volume fraction. The typical domain has side length L = 30-120 l in 2D and L = 10-30 l in 3D, where l = 0.25 in 2D and l = 0.35 in 3D, reducing any finite size effects. l is the non-dimensional length of a particle chosen so that the isolated swimming speed would be V 0 = 1. Due to the semidilute assumption, the flow acting on a particle is the sum of flows from all other particles including periodic images outside the computational domain. The periodic boundary conditions allowing for these flows are taken into account by generating virtual image particles and computing the flow due to each up to a distance from the reference cell beyond which the results do not change.
The present model captures the far-field hydrodynamic interactions while also introducing a mechanism for the displacement and realignment due to collisions. In addition, tumbling was introduced (see section A.1 of the appendix) for a more detailed comparison between experiment and simulations. Numerically, tumbling is accounted for in the two simulations depicted in figures 5 and 6 by choosing a random portion χ of the total number of particles to tumble in a given time interval t. In 2D, the rotational diffusion is then D θ = χπ 2 /3 t and in 3D D θ = 2χ π 2 /3 t (variance of the uniform distribution in the interval of [−π, π] is π 2 /3). Additional more complex steric effects such as lubrication forces may also contribute to the flow, but add greatly to the computational complexity and will not allow for the efficient simulation of a large number of particles.
In all cases, the initial orientation of each particle is taken to be random and the initial positions are uniformly distributed throughout the domain, except the situation considered in figures 2(a)-(d) where the bacteria were aligned along a large vortex. The simulations are then run until t = 40t = 10 s before any physical quantities are computed to allow the particles sufficient time to interact. Heret = l/V 0 is the characteristic time scale representing the time for a particle to swim its length. Figure 4 illustrates that beyond t = 10 s the system has settled and no major changes in the collective dynamics are observed.

Results
In this section we provide results for both the quasi-2D suspension and the 3D model. First we compare the results of the quasi-2D simulations to the experiments of Sokolov et al performed on a thin film suspension of Bacillus subtilis [2]. Figure 5 depicts the results of simulations showing the relationship between the correlation length L and the volume fraction . There are two clear regimes: individual motion ( < 0.2) and collective motion ( > 0.4) connected by a steep transition consistent with experimental observation. In addition, we calculated the mean particle velocity In dilute suspensions where 1 the mean particle velocity is essentially the isolated swimming speed V ≈ V 0 , but as the concentration becomes large the simulations show that the mean particle velocity increases, due to hydrodynamic entrainment, by a factor of almost six with concentration consistent with experimental observation [2]. Figure 5 illustrates the effects of tumbling. The first observation is that the larger the tumbling rate the higher the critical concentration is for the transition to collective motion. This is due to the fact that tumbling disrupts the alignment of particles (which is a characteristic of the collective state): more particles are needed to increase the contribution from collisions (that increase alignment) in order to dominate the effect of tumbling. Second, even in the absence of tumbling, there exists a critical concentration for the onset of collective motion. This behavior can be interpreted as a manifestation of 'intrinsic noise' generated due to interactions between multiple bacteria. This is consistent with the observation in our previous work [18], that hydrodynamic interactions in a deterministic system of many interacting particles leads to the same effect as random tumbling. Also, it can be seen from figure 5 that the correlation length, and, consequently, the scale of the collective motion is independent of the domain size for L > 30.
In addition, our simulations show that in three dimensions the correlation length of collective motion is essentially constant as the volume fraction changes above a certain critical value crit ≈ 0.2-0.4, figure 6(d), in agreement with experiment as well as being qualitatively and quantitatively similar to the results of the quasi-2D simulations in figure 5. This value is below the prediction of [11] that the transition to collective motion occurs when the volume fraction defined by the hydrodynamic radius of the particle becomes of the order one. Figure 6(a) shows the effect of swimming speed V 0 on the correlation length for a quasi-2D suspension. We observe that in a suspension where particles do not tumble the correlation length remains essentially constant as the swimming speed increases, showing the correlation length is independent of V 0 . If particles are allowed to tumble, then at low swimming speeds the correlation length rapidly decreases to the size of a particle (individual random motion). For higher swimming speeds the effect of tumbling on the collective motion is less pronounced. This shows directly that tumbling is the main cause of the decrease in the correlation length rather than the slow-down of the particles due to tumbling 6 . The effect of tumbling was also examined in full 3D simulations and the results are given in figure 6(b). We find that when the tumbling rate is increased, we recover the qualitative behavior observed in [22] and in 2D simulations: a decrease in the correlation length at low swimming speeds.
We now consider the effect of the system parameters on the correlation time. In the simulations, we normalize the correlation time by the characteristic time scale t 0 = l/V 0 in place of the mean free timet 0 (t 0 = 1/V 0 cl in 2D andt 0 = 1/V 0 cl 2 in 3D) between collisions used in experiment, where c is the concentration [22]. In addition to being independent of the concentration, the inset to figure 5 demonstrates that the rescaled correlation time, T /t 0 , is independent of the swimming speed in the collective regime. This is consistent with experiment where the correlation time is rescaled by the mean free time and was observed to be constant as the swimming speed increased [22,23]. One limitation of the present model is the introduction of a third time scale, collision time, in addition to the characteristic time scale and the time scale for collective motion: the collisions are soft and not instantaneous, and, therefore, time elapses as particles collide and reorient. The collision time, t col , is determined by the stiffness of the Lennard-Jones potential near the equilibrium distance r 0 = 2 1/6 σ LJ and is of the order t col ≈ σ 2 LJ /(36 × 2 2/3 ε). A simple calculation provides an estimate for the collision time in terms of the characteristic time scale and the mean free time t 0 : t col ≈ 0.2t 0 -0.4t 0 . However, our study revealed that this effect does not change the qualitative behavior of the correlation time.
Now that we have established that the model captures most experimental observations, we want to investigate relationships which are difficult and have not yet been observed experimentally. One such example is the effect of the dipole moment, characterizing the strength of hydrodynamic interactions, which depends on the particle shape and aspect ratio. The relative dipole moment U 0 decreases for slimmer particle with the decrease of its crosssection. Thus, for a given swimming speed and concentration, by changing U 0 we can compare the relative contribution of hydrodynamic interactions to collisions. Figure 7 illustrates the dependence of the normalized correlation length (L/l) on the dipole moment U 0 for a fixed swimming speed and volume fraction in the collective state. In quasi-2D we see a very minor increase in the correlation length and in 3D the correlation length increases rapidly as the dipole strength decreases. These results indicate that the aspect ratio plays a critical role in the onset of collective motion in a suspension of microswimmers, consistent with the observations made in [12]. Another example is the effect of the aspect ratio on the collective dynamics of the suspension. This can be quantified by the Bretherton constant B = a 2 −1 a 2 +1 , where a is the aspect ratio of the major and minor axes of an ellipsoid. Simulations show that the Now we briefly compare our results with past works studying the collective state. First, figure 2 shows that the power spectrum from the experiment and simulations is maximized below k = 0.3 consistent with the prediction of the instability of the isotropic state at low wave numbers [10,11]. We find the transition to the collective state occurs at a volume fraction between = 0.2 and 0.4, which is consistent with the prior results of [14] whose analysis yields the prediction of a critical concentration of = 0.3. In addition, Hernández-Ortiz and co-workers [15,16] show an increase in the mean particle velocity by an order of magnitude similar to the results in figure 5, which predict a six-fold increase in the mean particle velocity matching more closely to experiment [2]. While these highlighted results are consistent with previous works, new results were introduced such as the consideration of the onset of collective motion as a function of time or the direct comparison of the effects of hydrodynamic interactions versus collisions by studying the effect of the dipole moment in figure 7.

Conclusions
In this work we have introduced a simple model for studying correlation properties of mesoscopic collective motion in bacterial suspensions. Numerical simulations of this model demonstrate that both long-range hydrodynamic interactions and short-range collisions between bacteria are crucial to obtain agreement with experiment. Specifically, we were able to capture the qualitative effects on the correlation length and correlation time of a bacterial suspension as the main system parameters are varied: swimming speed, concentration, tumbling rate and dipole moment. In addition, we presented new results from simulations and experiments on the time needed to establish the collective motion.
Theoretical modeling also allows for results which are difficult to obtain experimentally. Our model allows for the examination of the effect of particle shape/aspect ratio on the collective dynamics. The main result was that the collective state is independent of the concentration for c > c cr and the swimming speed. Therefore, the correlation length and time in the collective state are essentially controlled by the particle size and shape, as suggested by experiments [12,22].
Our work reveals the delicate balance between various physical mechanisms determining collective motion in the suspensions of bacteria, as well as a critical role of the geometrical confinement and dimension on its correlation properties. Our experiments and simulations suggest that the organization of collective motion in bacterial suspensions is governed by two competing mechanisms: short-range collisions and long-range hydrodynamic interactions. While collisions tend to increase the correlation length of collective motion, hydrodynamic interactions result in its decrease. This is demonstrated, for example, by figure 7. The decrease of the hydrodynamic dipole strength in 3D results in a sharp increase of the correlation length. At the same time, in quasi-2D, the hydrodynamic interaction are screened due to effects of the confining walls and play a minor role on the correlation length. This observation is consistent with simulations in [12] where collective motion in two dimensions was studied essentially without long-range hydrodynamic forces.
A pure hydrodynamic origin of the collective motion in bacterial suspension was suggested in a number of works [7-9, 11, 28]. These works demonstrated that the hydrodynamic interactions alone are sufficient to predict the instability of the isotropic state for suspensions of pushers and persistent turbulent-like behavior beyond the critical concentration. However, the resulting spatial scale of collective motion was determined essentially by the size of the integration domain. Our model suggests that in suspensions studied in the experimental work [22] the main mechanism responsible for collective motion is in fact collisions, while the hydrodynamic interactions decrease the local order (e.g. correlation length). Also, in such suspensions the spatial scale of collective motion does not depend on the size of the domain. This implies that the collective states formed in systems described by pure hydrodynamic interactions, like in [7-9, 11, 28], may have a fundamentally different form of organization from that currently observed experimentally in [22].
While our model is computationally efficient and illuminates many interesting phenomena it does have some limitations. We obtained good qualitative agreement with experiment; however, the numerical value for the correlation length is 3-4 times smaller. The main reason for this discrepancy is the use of spherical particles (symmetric Lennard-Jones potential), with aspect ratio incorporated only in the parameter B. The use of true anisotropic particles is highly desirable, but ultimately leads to a drastic increase of the computation time. It was shown by experiment that hydrodynamic interactions and collisions are proportional to the swimming speed. In our model this remains mostly true, but some deviation is seen due to the third time scale that arises related to the collision time due to the use of a Lennard-Jones potential to model collisions. This time scale is only important at high swimming speeds, but did not inhibit the analysis of the qualitative effect of certain physical mechanisms on the properties of the collective dynamics of the system. by the NIH grant 1R01GM104978-01. The work of SDR was supported by the DOE grant DE-FG02-08ER25862.

A.1. Description of the model
The particle equations of motion are determined by the balance of forces and torques on each particleẋ where u is the flow defined in the model section depending on whether we consider a quasi-2D or full 3D model. This term couples the fluid Stokes equations and particle equations of motion. The translational motion (A.1) is a generalization of Stokes drag law f = 6π ηlu from the balance of forces and (A.2) is Jeffery's equation for the orientation of a passive ellipsoid [29] with additional contributions due to interactions and collisions. Modeling pure hydrodynamic interactions may cause the point dipole particles to overlap. This will result in an artificial divergence of velocities. Thus, we introduce a truncated purely repulsive Lennard-Jones type potential, W = 4ε[(σ LJ /r ) 12 − (σ LJ /r ) 6 ] + ε for r < 2 1/6 σ LJ and zero otherwise where ε = 0.01 in simulations. This enforces excluded volume constraints/soft collisions through a contribution to (A.1) by the repulsive force F = −∇W . σ LJ is a parameter representing the characteristic length l, r is the distance between two particles and ε is a parameter which determines the strength of repulsion. In contrast to [14][15][16] where the evolution equations only govern the positions of the beads in their dumbbell model, the particle orientations here evolve according to explicit equations of motion due to the torques acting on a given particle (Jeffery's equations, (A.2)). By enforcing Jeffery's equations we ensure that the point dipoles interact with the fluid as ellipsoids. The shape/aspect ratio is accounted for explicitly through the Bretherton constant, B (0/1 sphere/needle). This results in each point dipole bacterium having a size and shape.
The vorticity and rate of strain induced by hydrodynamic interactions are ω j ( is a realignment term due to collisions between particles. This term reinforces the following interaction between the bacteria: if particles approach one another near parallel (d i × d j small) or near perpendicular (d i · d j small), then there is only a small torque. Also the coefficient |F(x i − x j )| ensures that the torque after a short period of time will reorient the particles away from each other as seen in figure A.1. Here ξ controls the swim off rate and is taken to be 1 in the simulations presented in this work. This behavior is motivated by experimental observations of bacterial collisions [17] and the analysis of hydrodynamic interaction between two elongated self-propelled particles [30]. Simulations without this term yield the formation of densely packed clusters of bacteria separated by essentially empty domains. This behavior is in disagreement with experiment where a practically homogeneous distribution of bacteria is observed.  Since tumbling is present in the experimental data we also want to incorporate it into the simulations. Numerically the strength and frequency of tumbling, which is a random reorientation on the unit sphere, is controlled by picking a random percentage of the particles to tumble in a given predefined interval of time. From this D θ can be calculated for comparison with experiment. Unless otherwise noted all simulations are run in the absence of tumbling.

A.2. Comparison of quasi-two-dimensional and three-dimensional
In figure A.2 the streamlines associated with (a) quasi-2D and (b) full 3D models are depicted. One can see that they are qualitatively different and illustrate the effect of confinement on the suspension. Throughout this work the quasi-2D simulation results have been presented in detail and remarks are made indicating the differences with full 3D simulations when appropriate.

A.3. Correlation functions
We consider the spatial and time autocorrelation functions in the suspension by computing the coarse-grained flow field, V(x), generated by the fluid velocity, u(x), on a discrete grid. Once this field is computed the spatial correlation function can be defined as a function of distance r = |r| between points in the fluid where · x is an average over spatial coordinates. This is consistent with its calculation in experiment [22]. The correlation length is then extracted by fitting the time-averaged correlation function to an exponential function of the form e −r/L + A where A 1 and L is the correlation length as in [2]. Similarly the time correlation can be extracted from which is then spatially averaged. We non-dimensionalize the correlation length and time by the particle length l ∼ σ LJ and the characteristic time scale t 0 = l/V 0 (time for a bacterium to swim its length). Specifically, the first set of results in the results show the effect of the concentration and swimming speed on the correlation length and the latter results illustrate how the correlation time depends on these two parameters.