Simulations of Brownian tracer transport in squirmer suspensions

In addition to enabling movement towards environments with favourable living conditions, swimming by microorganisms has also been linked to enhanced mixing and improved nutrient uptake by their populations. Experimental studies have shown that Brownian tracer particles exhibit enhanced diffusion due to the swimmers, while theoretical models have linked this increase in diffusion to the flows generated by the swimming microorganisms, as well as collisions with the swimmers. In this study, we perform detailed simulations based on the force-coupling method and its recent extensions to the swimming and Brownian particles to examine tracer displacements and effective tracer diffusivity in squirmer suspensions. By isolating effects such as hydrodynamic or steric interactions, we provide physical insight into experimental measurements of the tracer displacement distribution. In addition, we extend results to the semi-dilute regime where the swimmer-swimmer interactions affect tracer transport and the effective tracer diffusivity no longer scales linearly with the swimmer volume fraction.


Introduction
Active suspensions of micro-swimmers such as spermatozoa (Creppy et al. (2015)), bacteria (Wensink et al. (2012)) or microalgae are common both in the natural environment, such as oceans (Pedley and Kessler (1992); Stocker (2012); Durham et al. (2013)), lakes, ponds and within living organisms, such as the human body. Such suspensions can exhibit the formation of coherent structures or complex flow patterns (Saintillan and Shelley (2011); Marchetti et al. (2013); Wensink et al. (2012)) which may lead to enhanced mixing of chemicals in the surrounding fluid, the alteration of suspension rheology (Rafaï et al. (2010); López et al. (2015)), or increased nutrient uptake. Mixing and transport of microscopic, inert 2 of 19 B. DELMOTTE ET AL.
particles by motile microorganisms have been a topic of recent interest as such suspensions are a prime example of out-of-equilibrium systems. The particles experience Brownian motion due to their small size and are further affected by hydrodynamics interactions and collisions with the micro-swimmers. Enhanced particle transport in active suspensions has been observed in the presence of collective motion (Wu and Libchaber (2000)), but also in the absence of it (Kurtuldu et al. (2011)). Understanding the mechanism of enhanced tracer transport can provide insight into biological processes such as predatorprey interactions and the plankton food chain (Kiørboe et al. (2014)), as well as chemical signalling and quorum-sensing (Kim et al. (2016)). In addition, the underlying mechanism may also provide the foundation for the design of novel biomimetic micro-fluidic devices that use similar strategies for enhanced mixing and stirring at small scales.
Generally speaking, the motion of a small particle in a suspension of micro-swimmers results from the interplay between Brownian diffusion due to thermal fluctuations and transport by the flow field induced by all the swimmers. The diffusion coefficient of the tracer, D 0 is set by the Stokes-Einstein diffusivity where ζ is the friction coefficient of the tracer particle and k B T the thermal energy. The tracer particles are advected by the velocity field u(x) which corresponds to the disturbances generated by the collection of micro-swimmers in the vicinity of the tracer. The velocity u(x) depends highly on the magnitude and the spatial rate of decay of the disturbances, as well as the local swimmer volume fraction, φ v . As the swimmers are moving, over time, the ever changing advection of the tracer particles results in their diffusive behaviour. In the dilute regime where φ v 1, many experiments (Leptos et al. (2009);Miño et al. (2011); Jepson et al. (2013); Kasyap et al. (2014)) have shown that the scaling between the tracer diffusivity, D act , due to swimmer activity and swimmer volume fraction is linear where D e f f is the total effective tracer diffusivity and α has the units of a diffusion coefficient. In terms of swimmer number density, n and swimming speed, U, this can also be expressed as D act = nUΛ , where nU is the "active flux" of micro-swimmers and Λ scales like the fourth power of the swimmer size based on dimensional analysis. In the dilute limit, where swimmers move along straight paths and the interactions between tracers and swimmers are well characterised, the effective diffusion coefficient can be computed (Lin et al. (2011);Miño et al. (2011);; ; Kasyap et al. (2014); Thiffeault (2014); Burkholder and Brady (2017)) by averaging the displacements of a single particle due to repeated interactions with swimmers that move independently of one another. In many of these studies, the swimmers are modelled as spherical squirmers. The spherical squirmer model was first introduced by Lighthill (1952) and furthered by Blake (1971) and provides the motion and flow generated by a spherical swimmer propelled by small distortions of its surface. Theoretical studies have shown that D act can be separated into two contributions: random swimmer reorientations and the entrainment by the swimmer along the straight trajectories. In the experiments of Leptos et al. (2009), the tracer diffusion coefficient in 3D dilute suspensions of C. Rheinardtii was obtained from data taken over period of two seconds. Since the reorientation time of C. Rheinardtii due to phase slips between flagella is approximately ten seconds ), swimmer reorientation does not significantly affect the tracer diffusion. Hence, by considering entrainment only,  obtained an estimate ((D e f f − D 0 )/φ v 83µm 2 s −1 ), surprisingly close to the experiments of Leptos et al. (2009) ((D e f f − D 0 )/φ v 81.3µm 2 s −1 ). Using a similar approach, Thiffeault (2014) 3 of 19 showed that the distribution of tracer displacements from Leptos et al. (2009) can also be reproduced using the squirmer model. However, these two works only consider point tracers while particles in experiments are micron-sized beads (Wu and Libchaber (2000); Leptos et al. (2009);Miño et al. (2011);Kasyap et al. (2014)), macromolecules, or dead cells (Jepson et al. (2013)) with a finite size. Recent experiments (Jeanneret et al. (2016)) in the dilute regime have shown that the front-mounted flagella of C. Rheinardtii can trap such micron-sized particles and generate very large displacements through direct entrainment. Such dramatic events are due to near-field interactions and are strongly related to swimmer geometry, as well as actuation mechanism.
Beyond the dilute limit, there can be a break down in the linear dependence of the effective diffusion coefficient on swimmer volume fraction. The departure from linearity, however, appears to depend strongly on the details of the system used for study. For example, Wu and Libchaber (2000) showed that the linear scaling holds up to φ v = 10% in two-dimensional films of E. coli, while Kasyap et al. (2014) observed in a three-dimensional bath that the linear trend breaks down for φ v 2.5%. As stated by Kasyap et al. (2014), the reason for the deviation from linearity is still unclear, and "one possibility is the occurrence of multi-bacterial effects on the tracer diffusivity." While suspensions of bacteria are known to display large-scale collective dynamics (Wu and Libchaber (2000); Dombrowski et al. (2004)) characterised by swirls and jets, the nonlinear dependence is also observe in swimmer suspensions that do not exhibit larger scale motion. For example, Kurtuldu et al. (2011) measured tracer displacements in a suspension of C. Rheinardtii confined to a thin film of liquid and obtained a power-law for the effective diffusivity D e f f /D 0 ∼ φ 3/2 v . A numerical investigation of the semi-dilute regime using the squirmer model and Stokesian Dynamics Ishikawa et al. (2010)) to compute the motion of non-Brownian fluid particles yielded a linear scaling of the tracer diffusivity with volume fraction for values up to φ v = 15%.
In this paper, we present results from simulations exploring tracer transport in dilute and semi-dilute suspensions of squirmers. In our simulations, we employ recently developed numerical tools based on the force-coupling method (FCM), (Maxey and Patel (2001); Lomholt and Maxey (2003); Keaveny (2014); ; ) that allow for multibody hydrodynamic interactions between active and passive particles, polydispersity in particle size, particle Brownian motion that satisfies the fluctuation-dissipation theorem, and steric interactions. A description of the squirmer model and fluctuating FCM are provided in Section 2. In the dilute regime, we obtain quantitative agreement between our simulation results and the experimental results of Leptos et al. (2009) for the effective tracer diffusion coefficient, as well as for the tracer displacement distribution. By selectively removing the flow disturbances due to the squirming modes, we quantify the contributions of hydrodynamic and steric interactions on tracer displacement and provide insight into physical mechanisms giving rise to particular features of the tracer displacement distribution. These results are shown in Section 3.1. We extend these results to semi-dilute concentrations in Section 3.2. Here, we examine the non-linear dependence of the effective tracer diffusion coefficient on the swimmer volume fraction, as well as characterise the distribution of tracers in the suspension. Finally, we discuss our results and future directions in Section 4.

Mathematical model for the simulations
In our simulations, we consider N p squirmers dispersed in fluid containing N t tracers giving a total number of particles N = N p + N t . All particles are spherical, however the squirmers have radius a sw , while the smaller tracers have radius a. The position of particle n is denoted by Y n , while its orientation is p n . The motion of both the spherical tracers and squirmers is governed by overdamped Langevin dynamics, which may be expressed as where Y is the 3N × 1 vector containing all particle positions, i.e.
. . , p T N ] T is that for the particle orientations. The matrix Q is the block diagonal matrix whose non-zero entries are given by Q 3n+k,3n+l = ε klm p n m for n = 1, . . . , N. The vectors V and W are the translational and angular velocities, respectively, of the particles that are related to the vector of forces, F , and torques, T , on the particles through where M is the 6N × 6N low Reynolds number mobility matrix for all particles. In Eq. (2.3), we indicate explicitly the four 3N × 3N submatrices relating either forces or torques with translational or angular velocities. In addition to V and W , the particles have velocities V sq and angular velocities W sq which encapsulate the swimming velocity, Up n , of each squirmer, as well as the additional velocities and angular velocities of all particles due to the flows induced by each of the first two squirming modes (Blake (1971); Ishikawa et al. (2006)). For a squirmer centred at the origin and in a frame moving with the squirmer, this induced flow is given by where B 1 is related to the swimming speed through B 1 = 3U/2 (Blake (1971)) and B 2 controls the strength and sign of the swimming stresslet, G, through The squirming parameter β = B 2 /B 1 describes the relative stresslet strength (Ishikawa et al. (2006)). For β > 0, the squirmer is a 'puller,' bringing fluid in along p and expelling it laterally, whereas if β < 0, the squirmer is a 'pusher,' expelling fluid along p and bringing it in laterally. The random velocities,Ṽ and angular velocities,W , obey the fluctuation-dissipation theorem and their statistics are related to the mobility matrix M through Also appearing in the equations of motion are the thermal drift terms k B T ∇ Y · M V F , k B T Q∇ Y · M W F , and −2k B T D W T P that arise from taking the overdamped limit and are required to obtain particle distribution dynamics that are governed by 2.1 Force-coupling method To compute the particle dynamics, we rely on the force-coupling method and two recent extensions to active particles ) based on the steady squirmer model (Lighthill (1952); Blake (1971); Ishikawa et al. (2006)) and to Brownian suspensions (Keaveny (2014); ) using concepts from fluctuating hydrodynamics. The force-coupling method uses regularized force distribution and volume averaging to generate the far-field approximation of the mobility matrix, M .
To compute the particle velocities and angular velocities arising from the forces and torques on the particles, as well as those due to squirming, we first solve for the Stokes flow where F n is the force particle n exerts on the fluid and S n is the stresslet of particle n due to its rigidity. We note in our simulations that all particles are torque-free, i.e. τ n = 0 for all n, and we ignore the tracer stresslets due to their small size. The only forces we consider are short-range pairwise repulsive forces that represent contact forces between two particles, and prevent overlap ). Additionally, for the squirmers, we have as to yield flow corresponding to Eq. 2.4 for each squimer. In Eq. (2.9), we also have the two Gaussian envelopes that are used to project the particle forces onto the fluid, The length scales σ n;∆ and σ n;Θ are related to the radius, a n , of particle n through σ n;∆ = a n / √ π and σ n;Θ = a n / 6 √ π 1/3 . After solving Eq. (2.9) for u, the velocity of each tracer is determined from while for the squirmers, the velocities, angular velocities, and local-rates-of-strain are given by where the terms W n and K n are included to subtract away artificial self-induced velocities and localrates-of-strain due to the volume integration of the squirming modes. Expressions for these terms are provided in . The local rate-of-strain for each squirmer is required to be zero, E n = 0, giving rise to the stresslets, S n . To compute the random velocities and angular velocities, we consider the Stokes flow where P is a fluctuating stress driving the random fluid flow. The statistics for P, in index notation, are given by (2.20) Using the FCM volume averaging operators, we first enforcẽ to obtain the values ofS n and subsequently computẽ As demonstrated in Keaveny (2014), the resulting velocities and angular velocities will satisfy the fluctuation-dissipation theorem with the covariance given by the FCM approximation of the mobility matrix. We note that while we have described the computation of deterministic and stochastic velocities separately, due to the linearity of Stokes equations, these can be combined into a single computation of the total particle velocities. While this describes how we obtain the deterministic and random motions of the particles, it remains to integrate the equations of motion with the correct thermal drift. We accomplish this by employing the midpoint drifter-corrector (DC) time integration scheme ) that inherently accounts for the drift terms without having to compute them directly. In words, the DC recovers the drift by initially advancing the particle positions to the midstep using only the random velocities determined from the fluctuating fluid flow with no stresslets for particle rigidity. At the midstep, the full computation is performed to find the deterministic velocities, as well as new random particle velocities using the same initial realisation of the fluctuating stress. The result is correlated motion between the initial and midstep velocities that then can be exploited to capture the correct thermal drift. For details of the scheme, including error expansions showing explicitly that the drift is recovered, the reader is referred to .

Simulations parameters and set-up
In our simulations, we set our parameters to match the C. rheindardtii experiments of Leptos et al. (2009). We set the swimmer to tracer radius ratio to a sw /a = 5 and set the swimming speed, and hence 7 of 19 B 1 , to match U = 100µms −1 = 20a sw s −1 . Following Thiffeault (2014), we set β = 0.5. The simulations are carried out in a triply periodic domain with edge length L = 23a sw . The Stokes equations are solved using a Fourier spectral method with N g = 192 grid points in each direction. The number of tracers is set to N t = 1255, resulting in a tracer volume fraction of less than 0.34%. The number of swimmers is varied from N p = 12 to 451 to examine very dilute swimmer suspensions of φ v = 0.4%, as well as semi-dilute cases where φ v = 15%. Finally, the thermal energy k B T is set to obtain the same distribution of tracers displacements as Leptos et al. (2009) in the absence of swimmers. We note that the resulting value of the diffusion, D 0 = 0.34µ 2 s −1 is slightly higher than the value D 0 = 0.28µ 2 s −1 given in Leptos et al. (2009). A snapshot taken from a representative simulation is shown in Figure 1.

Dilute regime
3.1.1 Tracer displacements In their experiments, Leptos et al. (2009) measured the time-dependent probability distribution function, P(∆ x, ∆t), for tracer displacements for suspensions with swimmer volume fractions φ v = 0−2.2% over 0.3s. As the beat period for the flagella of C. rheindardtii is T = 0.02s, this observation time corresponds to 15 beat cycles. They find that unlike previous experiments using bacterial baths (Wu and Libchaber (2000)), the tracer displacement distributions exhibit non-Gaussian tails. The non-Gaussian tails are attributed to rare entrainment events that occur when a tracer particle comes in close proximity to a swimmer's surface. We performed simulations corresponding to the same swimmer volume fractions and observation times as in the experiments of Leptos et al. (2009). The resulting tracer displacement distributions from our simulations are shown in Fig. 2a. Our squirmer simulations adequately capture the Gaussian core of the distributions. As in Leptos et al. (2009), we find that the distributions are self-similar with respect to the diffusive scaling ∆ x/(2D e f f ∆t) 1/2 (Fig. 2b). We find, however, that the larger, but rarer, displacement events related to the tails are slightly underestimated by the model. This is consistent with similar findings of tracer displacements by squirmers (Thiffeault (2014)). We attribute this to the fact that in the near-field, the squirmer does not replicate the flow induced by swimming C. Rheinardtii and the tails of the PDF for short-time tracer displacements depend on the details of the flow near the swimming micro-organism. The differences in entrainment due to differences in swimming behaviour have been examined previously in . Fig. 3 shows the distribution tails in more detail. Though over less than a decade of displacements, we observe that our simulation data follows a x −4 power-law, while the data from Leptos et al. (2009) appears to behave as x −3 . We note that Leptos et al. (2009) originally fit the non-Gaussian tails by exponentials. The x −4 decay was observed in experiments by Kurtuldu et al. (2011) and predicted by Pushkin and Yeomans (2014) for suspensions of dipolar swimmmers with random orientations.
Previous studies (Zaid et al. (2011);Thiffeault (2014)) have argued that the non-Gaussian tails are due to the shortness of the observation time coupled with the rarity of the entrainment events. To verify this assertion, we have computed the evolution of the tracer displacement PDF over a longer time ∆t = 0 − 2s as shown in Figure 4a for φ v = 2.2%. As ∆t increases, we see that the Gaussian core of the distribution broadens. We also observe that small shifts in the mean value appears, which we attribute to statistical fluctuations in the data. After rescaling by the diffusive timescale, we see in Figure 4b that the tails decay more rapidly and the distribution approaches a Gaussian as ∆t increases. This convergence to a Gaussian is observed for any volume fraction φ v = 0 − 2.2%, see Figure 4c, with the rate of convergence increasing as φ v increases. Thiffeault (2014) derived the criterion for reaching a Gaussian distribution. He showed that for spherical squirmers with β = 0.5, the time to      : reach Gaussianity is ∆t = 3.57, 1.74, 0.8, 0.5s for φ v = 0.4, 0.8, 1.6, 2.2% respectively. Our simulation data (Figure 4) are in inaccordance with these results. We also note that convergence to a Gaussian distribution in dilute suspensions was observed in experiments where the algal suspension was confined to a liquid film (Kurtuldu et al. (2011)). While our results point towards convergence to a Gaussian distribution, we note that in Zaid et al. (2011), using point sources with no steric interactions, they found that Gaussianity should arise only when the disturbance flow induced by the swimmers decays as r −n with n = 1, and if n 2, as is the case of squirmers, the fluid velocity distribution, and thus tracer displacements, should deviate from Gaussianity, even at long times. Their results are supported by the work of Rushkin et al. (2010) on dilute suspensions of Volvox (φ v 1.5%). Their study shows that when considering only settling forces (n = 1), the fluid velocity fluctuations follow a normal distribution. When accounting for the degenerate quadrupole due to ciliary beating (n = 3), the fluid velocity fluctuations exhibit strong deviations from Gaussianity, as confirmed by experimental data.

Enhanced Diffusion
In addition to the tracer displacement distribution, Leptos et al. (2009) measured the mean-squared displacement of tracers as a function of time and observed that the motion of tracers is diffusive obeying ∆ x 2 = 2D e f f t for all times. Figure 5 shows the mean-squared tracer displacement from our squirmer simulations. As φ v increases, we observe that the onset of the diffusive regime occurs after a period of anomalous transport at short times where ∆ x 2 ∼ t Θ , 1 < Θ < 2. We note that similar behaviour was observed in other experiments (Wu and Libchaber (2000); Kurtuldu et al. (2011); Kurihara et al. (2017)) and the theoretical prediction Θ = 3/2 has also been proposed (Kurihara et al. (2017)). Once we have reached the diffusive regime at ∆t = 2s, however, our values for the mean-squared displacement are very close to those given in Leptos et al. (2009). We extract the effective diffusion coefficient from our data by fitting the mean-squared displacement by the solution of the Langevin equation (Wu and Libchaber (2000)) where τ is the timescale over which ballistic motion ( ∆ x 2 ∼ 2 D e f f τ t 2 ) transitions to diffusive behavior ( ∆ x 2 ∼ 2D e f f t). In Figure 6, we compare the effective diffusion coefficient from our simulations with those obtained by Leptos et al. (2009). We can see that the simulated values match the experimental ones within the statistical errors reported for the experiments.
3.1.3 Steric vs. hydrodynamic interactions In our simulations, and in all experiments, the tracer particles have a finite size and experience contact forces with nearby swimmers. Hence, three phenomena contribute to the tracer diffusivity in the dilute regime: the flows generated by the swimmers, collisions with swimmers due to contact forces, and tracer Brownian motion. In order to quantify the effect of each, we consider three situations: (i) a bidisperse suspension of passive particles achieved by setting U = 0, H n = 0, and G n = 0 for all swimmers, (ii) a suspension where the swimmers move in the fluid without generating any swimming disturbances and interacting with the tracers through contact forces (U = 0, H n = 0, G n = 0), and (iii) the full simulation model. Figure 7 shows the resulting PDF of tracer displacements for a dilute suspension (φ v = 2.2%) of squirmers with β = 0.5. In the purely passive case (i), the PDF is Gaussian as expected. For case  φ  (ii), the motion of the larger particles and collisions with the tracers generate a non-Gaussian PDF. The Gaussian core of this PDF is similar to that observed in case (i), but we see that there are additional non-Gaussian tails with power-law decay. When moving to the full model (case (iii)), the Gaussian core widens and the power-law decay of the tails increases. This widening of the Gaussian core was also seen in the experiments. These results show that rare swimmer-tracer collisions produce the large displacements that lead to the non-Gaussian tails, while the flows generated by the swimmers act to broaden the Gaussian core, thus enhancing the diffusivity of the tracer particles at short times.
3.1.4 Role of time-dependence The flows generated by C. rheinardtii are time-dependent (Guasto et al. (2010)) due to the beating of its two flagella. We can include this time dependence by allowing the parameters B 1 and B 2 appearing in the steady squimer model to be periodic functions of time. These functions can then be tuned to match the time-dependent swimming speed of C. rheinardtii and the location of the flow stagnation point as measured in Guasto et al. (2010). The details of our tuning procedure can be found in . After tuning, we find the average value of the squirming parameter over one beat period isβ = 0.12. Figure 8 shows the PDF of tracer displacements for φ v = 2.2% from the fully time-dependent model. For comparison, we have also included results from simulations of the time-dependent model with β = 0.5, as well as results from the steady model with β = 0.12 and β = 0.5. For β = 0.12, we find that the steady and time-dependent PDFs are nearly indistinguishable. Asβ increases toβ = 0.5, the PDF is slightly narrower with respect to that of its corresponding steady simulation β = 0.5. The difference, however, is quite small, and we find that time-dependence of the squirming modes does not significantly affect the PDF of the tracer displacements. These results are in accordance with Dunkel et al. (2010) that :β = 0.5, time-dependent. : β = 0.5, steady.
showed in the dilute regime the time-dependent and stroke-averaged swimming disturbances provide similar tracer scattering.

Semi-dilute regime
In this section, we perform simulations on more concentrated suspensions using the same value for the squirming parameter β = 0.5. Figure 9 shows the time-dependent PDF for tracer displacements P(∆ x/a, ∆t) at times ∆t = 0.3s and ∆t = 2s, and for volume fractions φ v = 0 − 15%. For all volume fractions, we observe that the PDF tends towards a Gaussian at long times, however, the convergence rate increases with volume fraction. Using an analysis based on singular flows fields, Zaid et al. (2011) state that the tracer displacement distribution should only be Gaussian for φ v > 25% if the flows generated by the swimmers decay like r −n with n 2. The appearance of a Gaussian, however, agrees with the theoretical predictions of Thiffeault (2014), as well as the experimental observations of Kurtuldu et al. (2011) that observed an increase in Gaussianity with concentration. This, perhaps, highlights the importance of resolving the finite-size of both swimmers and tracers to the observation of a Gaussian distribution of displacements. In addition, we note a shift of the mean at long times for high volume fraction. This is attributed to the onset of polar ordering ) that is typically observed in periodic squirmer suspensions.
As in the experiments of Kurtuldu et al. (2011) and Kasyap et al. (2014), our simulations show ( Figure 10a) the linear scaling (D e f f − D 0 )/φ v = const. breaks down for φ v > 2.2%, though we do not observe a clear power-law as in Kurtuldu et al. (2011). Additionally, our simulations yield D e f f /D 0 = 11 − 22 for φ v = 5 − 10%, while Kurtuldu et al. (2011) measured D e f f /D 0 900 for φ v = 7%. The large difference in enhanced diffusion can be due to our simulations being three-dimensional, while Kurtuldu et al. (2011) carried out experiments into a two-dimensional film with thickness on the order of the FIG. 9: PDF for tracer displacements at times ∆t = 0.3s (a), and ∆t = 2s (b), for φ v = 0 − 15%. cells' diameter (H 15 ± 5µm), resulting in a slower decay of the flows induced by the swimming cells, as well as an increased frequency of swimmer-tracer collisions. To quantify the effect of swimmer concentration on the swimmers' motion, we examine the swimmer mean-squared displacement as a function of concentration, as shown in Figure 10b. While we see a slight decrease in the slope as φ v increases, we observe that the swimmers are ballistic for all concentrations. Longer simulation times would be required to reach a diffusive regime in order to determine swimmer effective diffusivity D e f f ,sw , as done in Ishikawa et al. (2010). Our simulations, therefore, correspond to a different regime than that explored in Ishikawa et al. (2010) where both swimmers and tracers were diffusive. We note, however, that the values |β | > 1 used in Ishikawa et al. (2010) are greater than the values used here to match experimental data and higher values of β do lead to the onset of diffusive behaviour for the swimmers at earlier times.
While for each concentration the swimmers' motion is ballistic, we do observe clear differences in the tracer distribution as the swimmer concentration increases. Figure 11 shows the squirmer-tracer pair distribution function g(r, θ ) for φ v = 2.2 and 15%. The distribution function g(r, θ ) corresponds to the likelihood, relative to a uniform distribution, of finding a tracer at a distance r from the swimmer centre of mass and with its relative position vector forming an angle θ with the swimmer orientation. At higher concentrations, there is a much higher probability of finding a tracer directly in front of the swimmer, while the probability of finding tracers aft of the swimmer is reduced. This suggests that accurately representing the details of swimmer-tracer interactions for tracers directly in front of the swimmer is important in reproducing the large displacements. In addition, the increased localisation of tracer particles in the vicinity of the swimmers that occurs at higher swimmer concentrations could lead to increased rates of nutrient uptake by the population.

Discussion and conclusion
In this study, we performed simulations of tracer motion in dilute and semidilute suspensions of swimming particles. Our model includes the coupled effects of the finite tracer size relative to that of the swimmers, particle Brownian motion satisfying the fluctuation-dissipation theorem, and the flow disturbances induced by the swimmers using the steady squirmer model. For dilute suspensions, our simulations reproduce quantitatively the tracer displacement distribution and tracer diffusivity measured experimentally by Leptos et al. (2009), as well as those predicted theoretically by Thiffeault (2014). We demonstrate that the non-Gaussian tails of the tracer displacement distribution are linked to collisions between the swimmers and tracers, while the width of the Gaussian core is set by the many-body hydrodynamic interactions between the tracers and swimmers, as well as latent tracer diffusion. We also show that time-dependence of the squirmer modes has little effect on the distribution. In addition, our results demonstrate that at longer observation times, the number of swimmer-tracer collisions increase leading to the disappearance of the non-Gaussian power law tails and a Gaussian tracer displacement distribution. This corresponds with the predictions of Thiffeault (2014). In the semi-dilute regime, the swimmer-tracer collision rate increases and we observe a faster convergence to a Gaussian distribution of displacements. Additionally, the simulations yield a nonlinear dependence of the effective diffusion coefficient on φ v though the enhancement observed was much more modest than that measured (Kurtuldu et al. (2011)) in the suspensions confined to a thin film. In the nonlinear regime, we found that at the simulation timescales, the swimmers were still behaving ballistically, but there was a notable increase in the likelihood of finding tracers directly in front of the swimmers.
A particularly interesting property of this system is the importance of near contact swimmer-tracer collisions, even in dilute systems, as they lead to rare, but very large displacements. Even though the squirmer model provides a reasonable description of the far-field flow of many microorganisms, it is the near-field details that matter when resolving the collisions ). Thus, at longer times when there are sufficiently many collision events, the effective tracer diffusion coefficient will be determined almost exclusively by the collisions. The details of how the organisms generate flow and possibly deform their surfaces or move their flagella are then crucial in quantifying the effective diffusion (Jeanneret et al. (2016)). This could have important implications in predator-prey interactions, making only certain locomotion strategies viable for particle capture. For example, a swimmer with a rigid front surface could possibly be more effective at carrying along small particles as it swims, allowing more time for them to be ingested. In addition, the importance of near-field interactions on particle diffusion may certainly impact the design of synthetic artificial swimmers or other colloidal active particles for tracer mixing in microfluidic devices. Exploring the impact of different near-field swimmer-tracer interactions for swimmers that induce the same far-field flow, as well as how heterogeneity throughout the swimmer population affect tracer transport are areas of interest for future investigation.
In bacterial suspensions, the nonlinear dependence of the effective tracer diffusion coefficient on swimmer volume fraction has been attributed (Kasyap et al. (2014)) to the onset of large-scale collective motion in the suspension. In squirmer suspensions, especially at semi-dilute concentrations, the longtime dispersion properties might well be affected by polar ordering ), which may lead to net tracer transport along a particular axis, possibly aligned with the mean swimming direction.
Finally, it would be interesting to understand the impact of confinement on tracer diffusivity, especially with regard to the thin film experiments of Kurtuldu et al. (2011). Due to the stress-free boundaries, the hydrodynamic interactions are longer ranged than in bulk and may therefore enhance tracer diffusivity more than in unconfined suspensions. Additionally, we expect the swimmer-tracer collision rate to increase in confined systems which should further increase long time tracer diffusivity. Our numer-ical framework can simulate thin film geometries with stress-free boundary conditions while retaining Brownian motion and hydrodynamic interactions ). These simulations may provide some insight into the role of confinement and boundary condtions on tracer transport and form the basis of future studies.