Influence of thermal fluctuations on active diffusion at large P\'{e}clet numbers

Wavelet Monte Carlo dynamics simulations are used to study the dynamics of passive particles in the presence of microswimmers, taking account of the often-omitted thermal motion alongside the hydrodynamic flows generated by the swimmers. Although the P\'{e}clet numbers considered are large, we find the thermal motion to have a significant effect on the dynamics of our passive particles, and can be included as a decorrelation factor in the velocity autocorrelation with a decay time proportional to the P\'{e}clet number. Similar decorrelation factors come from swimmer rotations, e.g.~run and tumble motion, and apply to both entrainment and far field loop contributions. These decorrelation factors lead to active diffusivity having a weak apparent power law close to $Pe^{0.2}$ for small tracer-like particles at P\'{e}clet numbers appropriate for E. coli swimmers at room temperature. Meanwhile, the reduced hydrodynamic response of large particles to nearby forces has a corresponding reduction in active diffusivity in that regime. Together, they lead to a non-monotonic dependence of active diffusivity on particle size that can shed light on similar behaviour observed in experiments by Patteson et al.

Wavelet Monte Carlo dynamics simulations are used to study the dynamics of passive particles in the presence of microswimmers, taking account of the often-omitted thermal motion alongside the hydrodynamic flows generated by the swimmers. Although the Péclet numbers considered are large, we find the thermal motion to have a significant effect on the dynamics of our passive particles, and can be included as a decorrelation factor in the velocity autocorrelation with a decay time proportional to the Péclet number. Similar decorrelation factors come from swimmer rotations, e.g. run and tumble motion, and apply to both entrainment and far field loop contributions. These decorrelation factors lead to active diffusivity having a weak apparent power law close to Pe 0.2 for small tracer-like particles at Péclet numbers appropriate for E. coli swimmers at room temperature. Meanwhile, the reduced hydrodynamic response of large particles to nearby forces has a corresponding reduction in active diffusivity in that regime. Together, they lead to a non-monotonic dependence of active diffusivity on particle size that can shed light on similar behaviour observed in experiments by Patteson
While the precise results vary with the details of each system, all find the passive particles to exhibit enhanced motion over and above their own thermally-driven Brownian motion.
This swimmer-induced motion is super-diffusive on the time-scales of interactions with passing swimmers, and diffusive thereafter 7,14 .
Until recently the size of the passive particles has received little attention, with a range of sizes used across the literature but typically a constant size within a given study. Nevertheless one can identify several relevant properties that change with particle size: the thermal diffusivity and by extension the Péclet number Pe (defined as the ratio of advective and diffusive transport rates); the range of steric interactions with swimmers; and the near-field hydrodynamic response to swimmers [28][29][30][31] .
The influence of Péclet number was included in a theoretical study by Kasyap et al. 32 , which predicts a peak in swimmer-induced diffusivities at moderate Pe, rising from 0 at Pe = 0 (corresponding to infinitesimal tracers) and limiting to a finite value as Pe → ∞.
While this work accounted for particle size in the passive particles' thermal motion, they were hydrodynamically coupled to the swimmers by the (unregularised) Oseen tensor so that they were treated as point particles for hydrodynamic purposes. Their results might therefore not be expected to quantitatively match real systems except when at small Pe. Shum and Yeomans have performed detailed boundary element simulations of single swimmer-passive interactions with a wide range of passive particle sizes, neglecting all thermal motion 17 . From these results they obtained the dilute limit active diffusivity by integrating over impact parameters. In doing so they found a non-trivial dependence on the ratio of passive and active particle size, with a maximum active diffusivity for similarly sized particles when using a squirmer type swimmer model, but a minimum instead with bacteria models.
An experimental investigation of passive particle size active systems of E. coli was conducted by Patteson et al. 33 , who found a non-monotonic variation in active diffusivity, peaking when the swimmer and passive particle sizes are similar. This runs counter to the results for Shum and Yeomans' bacterial model. This non-monotonicity suggests regimes where different particle sizes are dominated by different physics. Since thermal diffusivity rises rapidly as particles get smaller, D ∼ R −1 , we hypothesise this is important and will study simulations of analogous systems with this in mind.
It is notable that thermal fluctuations have usually been omitted from simulations of microswimmers on the basis that swimmer Péclet numbers are generally much greater than unity. Work that has included thermal motion has done so in a hydrodynamically decoupled way 15 . While this can capture some of the physics, the hydrodynamic coupling is required to provide the correct relative thermal motion.
To elucidate the potential importance of thermal fluctuations we draw parallels to Taylor dispersion in pipes containing a steady shear flow 10,34 . In the absence of thermal fluctuations passive particles are carried parallel to the pipe axis at constant speed which is fastest at the centre of the pipe. Thermal fluctuations allow particles to cross stream lines, with a corresponding change of advective velocity, leading to dispersion of particles in the stream direction. Although the microswimmer flow fields are more complex, the crossing of stream lines is still expected to disperse advective motion, and will therefore have some influence on the active diffusion of passive particles, even if the Brownian motion itself is small. In the slightly more limited context of particle entrainment, this effect has been seen to lead to a non-monotonic distribution of entrainment jump sizes with particle size 10 .
Another reason for the limited microswimmer work including thermal fluctuations is the great computational expense required to include the correct hydrodynamics of thermal fluctuations, which is a challenge well known to the polymer community where it has spawned a wide range of simulation algorithms [35][36][37] . Here we make use of the recently developed Wavelet Monte Carlo dynamics (WMCD) algorithm to include hydrodynamically coupled thermal fluctuations efficiently 37,38 .
After setting out the hydrodynamic theory and simulation details in Sections II and III, we will demonstrate the validity of using WMCD for active systems in Section IV, where we use trajectories of passive particles in the flow of a swimmer on an infinite straight path as a test case. We then study the dynamics of passive particles in dilute mixtures of swimmers, looking in detail at the role of particle size and temperature on the velocity autocorrelation and active diffusivity in Sections V and VI respectively. In doing so we find we can include the effect of all relevant time scales through exponential decay factors in the velocity autocorrelation, allowing us to construct an ansatz function that successfully captures the complex behaviour seen in the active diffusivity.

A. Active diffusivity
The total diffusivity of a particle can be split into thermal and active contributions as The thermal diffusivity for a sphere of radius a in fluid of viscosity η is given by the wellknown Stokes-Einstein relation, which in a periodic cubic box of side length L is corrected to first order 39 .
Rather than thermal fluctuations, the active diffusivity D A is driven by hydrodynamic and steric interactions with active particles in the system, and can itself be written as the sum of those contributions 40 . To simplify data analysis we will not include steric interactions in this work so D A is purely hydrodynamic.
Regardless of its contributions, the active diffusivity can be expressed in terms of the velocity autocorrelation via the Green-Kubo relation 41,42 The lower limit 0 + denotes time t → 0 from above, such that the thermal contribution is excluded when working in the overdamped limit (see the next section). In practice this means the t = 0 value we use in our data is extrapolated back from data at small but finite t.

B. Equations of motion
Working on time scales where thermal fluctuations can be considered instantaneous, or equivalently on time scales longer than the fluid relaxation time, leads to the overdamped Langevin equations for translational (superscript T) and rotational (R) velocities: and where F j and Γ j denote the force and torque at r j , which may or may not be located on a particle. In this work i will correspond to a particle, while the sum over j corresponds to a sum over the swimmer forces described in Section III B. We will not apply any point torques, so Γ j = 0 for all j, but their inclusion here is useful for introducing the rotational mobility tensors which will be needed for correlations between the thermal fluctuations, ξ and Ξ.
The Lorentz reciprocal theorem links the TR and RT tensors by the transpose 43 For spheres, the unregularised versions of these tensors are to leading order in 1/r: where a i is the radius of particle i, δ ij is the Kronecker-delta and [r] × is the skew-symmetric tensor expressed as ε abcrb in index notation.
The mobility tensors used in this work are those that appear in Wavelet Monte Carlo Dynamics (WMCD) 37,38 , described below, which smoothly bridge the large r and δ ij terms.
Although reached in a completely different way, they closely approximate the tensors obtained by using Faxén's laws to sum the fluid flow over the particle surface 29 .
Finally, the fluctuation dissipation theorem gives 44,45 as the correlations between noise terms.

A. Wavelet Monte Carlo dynamics
We use a smart WMCD simulation algorithm, for which the full details can be found in The smooth approach to the r = 0 hydrodynamic tensors is achieved by setting finite values for minimum wavelet radii λ min , which is chosen separately for both translations and rotations at each particle size to give the appropriate particle mobility at r = 0.
The end result is an efficient algorithm that includes long ranged hydrodynamic correlations for both causal and thermal forces with a computational cost that rises with the total number of particles as N ln N per unit of physical time. The price is that it is limited to the mobility tensors for spheres and cannot currently handle lubrication forces or no-slip boundary conditions on the sphere surface. Nevertheless, the efficient inclusion of hydrodynamically coupled thermal fluctuations means WMCD is well placed to investigate whether these play a role in active-passive mixtures.

B. Swimmer model
Our swimmers are represented by a simple two-force pusher-type model, with a forward force F s placed at the swimmer centre and a tail force −F s placed A s = 3a s behind the swimmer, as depicted in Fig. 1. Using the WMCD mobility tensor these produce a swimming Encounters of passive and swimmer particles are governed by competition between the above swimming and the relative thermal diffusivity characterised by its large separation value We can define a corresponding Péclet number by using a s as the relevant length scale, and comparing the advection rate v s /a s to the diffusion rate D rel /a 2 s , leading to where Pe s ≈ 0.5F s a s /k B T is the intrinsic Péclet number of the swimmer alone.
The run and tumble motion exhibited by many micro-organisms, such as E. coli, is characterised by alternating phases of swimming and stopping, with increased rotational motion during the stopped 'tumble' phase. This has previously been modelled in simulations by instantaneous and random reorientations of swimmers at a Poisson distributed frequency 15,47 .
In this work the tumble phase is stretched out over a finite time t tumble , during which there is an increased rotational diffusion and the swimming forces are turned off, so that the swimmer temporarily becomes a passive particle. To reduce the number of variables and make the data easier to decipher, we use fixed values of t run and t tumble rather than choosing them from a Poisson distribution. Hence each swimmer in the system cycles between running and tumbling with the same period t run + t tumble , although each is initialised at a different point in this cycle. Table I lists the values of the physical parameters used for the results sections.

C. Simulation parameters
equalling 40 in our simulations. These runs dominate the time elapsed, so the active diffusion of the swimmers is well approximated by D

(s)
A ≈ a s v s λ/6. Because we have a finite t tumble we need to specify how much swimmer orientations decorrelate when tumbling. This is done via 48 which will be useful for relating tumble angles to decorrelation times in Section V. In our simulations we used θ tumble = arccos v(t tumble ) ·v(0) = 70 • , inside the range of angles identified in Ref. 47. Reorientation from tumbling and normal rotational diffusion are comparable when averaged over the run and tumble cycle, with tumbling being dominant and sub-dominant at the lower and higher temperatures respectively in Sections V and VI. We note that the small size of E. coli is important here, and thermal rotations would be less significant were we modelling a larger microswimmer.

The swimmer volume fraction in Sections
The volume fraction of passive particles is unimportant because they do not influence each other's motion, despite having correlated displacements. However, what is important for good statistics is the product of N p and total data collection time, which was a minimum of 1.85 × 10 5 t run of effective single particle tracking time per data point.

IV. LOW NOISE PARTICLE TRAJECTORIES
Our first results focus on interactions between individual swimmers and passive particles in the idealised scenario where swimmers move along an infinite straight path at constant speed. Such simulations have been done with more sophisticated techniques previously 16,17,19,40,49 , but are revisited here because they validate the use of WMCD for active systems while helping to visualise behaviour quantified in later sections.
Swimming along a perfectly straight path is not possible in WMCD, but a good approximation was achieved by switching off all rotations and reducing the temperature to raise the swimmer Péclet number to 1.3 × 10 8 . This also makes the Péclet number of relative motion high enough that thermal diffusion has negligible role in particle encounters.
In these simulations, performed in an infinite box, a single swimmer was set swimming in a straight line between −40a sx and +40a sx . A single passive particle was placed an impact parameter ρ off the swimmer's path at r p (t = 0) = ρŷ, and its position was traced out as the swimmer passed by. The x-y components of these trajectories are plotted in Fig. 2(a) for various ρ.
Qualitatively, these trajectories match expectations by forming (almost) closed loops with cusps 19 at large ρ, whilst at small ρ the cusps become more rounded and the loops open up with a significant finite net displacement. The details of how the loops change at small ρ are sensitive to the near-field details of the hydrodynamic mobility tensors and hence to passive particle radius a p , as evidenced by the clear differences between the loops for a p = a s and 2a s at same ρ. To demonstrate this sensitivity to particle size, Fig. 3 shows the flow fields experienced by passive particles of different sizes, as seen in the swimmer's reference frame. The most notable feature is the shaded recirculation zone (closed stream lines) close to the swimmer body that appears for small a p , which the passive particles do not enter as the swimmer passes by. The a p = 0.25a s tile is very close to the flow-field that would be produced using the unmodified Oseen tensor, while the differences in the other 3 tiles arise due to the near-field corrections in the WMCD tensor. These corrections are therefore responsible for the differences in the trajectories in Fig. 2(a). This highlights the limitations of treating passive particles as infinitesimal tracers when near-field flows are important. The trajectories for smaller ρ have positive net parallel displacements. For the limited but relevant range a p < ρ < A s we obtain theoretically the much weaker dependence ∆(ρ) ∼ +ρ −1 by treating the swimmer as two explicit point forces. In practice the crossover from positive to negative displacements can be seen in Fig. 2 This behaviour is akin to entrainment observed in real systems 9,10 , albeit driven by internal flows rather than steric repulsion and a non-slip boundary at the swimmer surface. We therefore describe such trajectories as 'entrained' in the following sections.

V. VELOCITY AUTOCORRELATIONS
Next we discuss the velocity autocorrelation of passive particles in active systems. The systems are now made periodic with cubic box, we use a reference temperature T 0 representative of 300K, and swimmers are free to rotate by rotational diffusion, hydrodynamic interactions and tumbling. Note that in these simulations we confirmed D A is proportional to swimmer concentration, as demonstrated in Fig. 4 whose results are discussed in detail in Section VI, so we are in the dilute limit where swimmer-swimmer interactions can be neglected. We also note that we find D A to vary roughly as (a p /a s ) 0.2 at small a p . This power cannot be explained by considering any one mechanism, and is a sign that we are in a complex regime with many contributing factors. D A itself has integrated out too much information to unpick these factors, motivating our focus on velocity autocorrelations instead.
We begin with a qualitative discussion of these and will use this as the basis for a quantita- tive discussion in Section V A. We will denote passive particle and swimmer autocorrelations with C vv we find in our data, alongside an example C (s) vv which is a simple exponential decay with some fine details coming from run and tumble motion which we will ignore.
In curve 1, which is typical of systems with a p a s , we see negative tails compatible with the forwards-backwards movement in the loop trajectories discussed in the preceding section. In curve 2, with smaller a p at the same temperature, this negative tail appears to have vanished, or at least been reduced to the size of noise in our data.
Note the value of C (p) vv (0) has increased between these curves, which is due to small particles having a larger response to the swimming forces in the near-field. The size of the increase, at more than a factor of 2, is indicative of how much the near-field contributes to C (p) vv (0), and we anticipate a similar rise in D A . Importantly, however, C (p) vv (0) reaches a maximum at around a p = 0.5a s , below which it is essentially constant. This will be discussed in greater detail in Section V A 2.
Finally, curve 3 in Fig. 5(a) shows that the reduction of temperature reveals a longtime exponential tail. This tail is present in the all our T = 0.1T 0 data with a p a s / √ 2, although its amplitude is not always the same. This tail runs almost parallel to C (s) vv in curve 4, suggesting the passive particles are tracking the swimmer motion, that is they are entrained. This is further supported by the quantitative similarity between C  There is in fact a subtle but important difference in the gradients of curves 3 and 4 at large times. First, tracing the long-time exponential back to t = 0 leads to a value larger than φ s C (s) vv (0), fitting with the entrainment volume around a swimmer being larger than the swimmer's own volume, as per Fig. 3. Second, this means the mid-time section of the curve undershoots the long-time exponential, which is consistent with there being a negative contribution akin to that seen in curve 1. This too is expected as the loop trajectories should still be present, and indeed the location of the peak of the negative (dashed) section of curve 1 coincides with the depression in curve 3. We believe this story applies to curve 2 as well, but in this case the negative contribution is closely matched to the positive tail so they cancel each other out.
Our final comment on the qualitative features of C (p) vv is that we find curves at different temperatures but the same particle size can be successfully mapped onto each other by introducing an exponential decay factor between them. This is demonstrated in Fig. 5(b), and will be central to our approach going forwards.

A. Ansatz function
We now attempt to quantify the velocity autocorrelations, starting with the swimmers since they have the simplest form, and we have seen that the passive particles can pick up the same form.
Ignoring some fine details seen in Fig. 5(a) that arise due to the fixed times in the run and tumble cycle, it is clear we have a simple exponential decay of the form, where the run and tumble factors account for the time spent not swimming and hence why the swimmer curve begins just below 1. The exponential decay comes only from reorientations because translational diffusion does not change the swimming direction, and therefore does not affect C (s) vv beyond the Brownian spike at t = 0, which we are ignoring. Hence we have where the decay time associated with normal rotational diffusion is and the decay time for run and tumble motion, spread across the whole run and tumble cycle, is Note the final expression uses Eq. (18) and our typical tumble angle of 70 • . Together, these predict τ s v s /a s = 39.2, in good agreement with the decay time observed in curve 4 in Fig. 5(a).
We also expect the C (p) vv to decay with τ s . In the entrainment tail the reasoning is the same as for the swimmers -the decorrelation of the direction of v(t) -but we also expect it to play a role in decorrelating non-entrained loops where swimmer rotations effectively force the passive particles onto different stream lines, as well as rotating the flow field. By similar reasoning, we expect an additional decorrelation time coming from (translational) Brownian motion, τ B . This should also feature in the entrainment tail where it drives entry and escape of the entrained volume. In higher density systems we might expect a similar term for swimmer-swimmer interactions, but we do not consider those here.
We express the total C (p) vv as the sum of terms from entrainment and non-entrained loops: which we now detail separately.

The entrainment term
Entrained particles follow swimmers, therefore we expect C (Ent) vv to be given by Eq. (19) modified to account for Brownian escape and the actual entrained volume fraction φ Ent .
This implies that The details of the decorrelation factor, including the perhaps unexpected square root, will be discussed after φ Ent .
As Fig. 3 showed, φ Ent varies with particle size, and has 3 regimes of behaviour. For a p > a s , there is no entrainment as the response to the swimming forces is always less than that of the swimmer. For a p ≪ a s the entrainment volume is constant and is approximately 4 times the swimmer volume for our model. As a p increases the volume begins to decrease when the hydrodynamic near-field of the particle, i.e the distance within which the hydrodynamic response is different to the Oseen tensor, is comparable to the geometric size of the swimmer, A s . In our WMCD simulations that hydrodynamic range is 5.35a p , leading us to define the ratio R = 5.35a p A s = 1.78 a p a s (25) with which we would anticipate the regime change at R ≈ 1. Assuming a simple linear interpolation between the large and small a p regimes, we have It is the asymmetry in this first passage problem that leads to the square root rather than the simple exponential decay used in other decay factors.
That calculation gives us a handle on the form of the decay time τ B,Ent ∼ a 2 s /D rel = (a s /v s )Pe rel . Using all the other parameters in C (Ent) vv as described above, we set the numerical prefactor to match the entrainment tail in Fig. 5(a) curve 3. This gives us corresponding to diffusing a distance of 6D rel τ B,Ent ≈ 0.65a s , which is reassuringly less than a s .

The loop term
The non-entrained loop contribution in Eq. (23) needs to provide the negative tail seen in the a p = 2a s data in Fig. 5. The functional form of this could in principle be calculated in the far-field where the loop trajectory is known mathematically for a dipole swimmer 19 , but the near-field is not entirely absent in our expression so we would expect model-dependent terms to enter. We therefore choose to take a simpler route and use a functional form that has the correct features: The task is now to identify all the new parameters, starting with the coefficients c 0 and c 2 . c 2 will not be explored in detail, but we note it must satisfy 0 < c 2 1/2 to ensure both a negative tail exists and C c 0 sets the value at t = 0 and can be written as The total C (p) vv (0) is most easy to access in our data by interpolating small t values back to 0 assuming a simple exponential decay, thereby avoiding the Brownian spike which overwhelms the zero-time data. Values obtained by this procedure are shown in Fig. 6(a), where we see vv is flat at small a p /a s and decays at larger values. The flat regime is simply the result of passive particles being small compared to the distance between the swimming forces, A s , so they act like infinitesimal tracer particles.
Indeed, the regime change occurs close to R = 1 ⇔ a p /a s = 0.56, supporting this picture.
What is less easy to understand is the apparent dependence on temperature in this regime, with the high temperature data being too small to be accounted for by our error margins.
We attribute this apparent T -dependence to the assumption of a simple exponential decay when interpolating back, which underestimates the contribution from exp[− t/τ B,Ent ]. The corresponding error is largest at small Pe rel and only when the entrainment term is present, both fitting with where the difference occurs in Fig. 6(a).
The decay at larger a p /a s can be understood using scaling arguments, detailed in Appendix B, which use the fact our mobility tensor can be written as a −1 p G(r/a p ) when a p ≫ A s , leading to c 0 ∼ a −1 p . We can capture both regimes with the piecewise function where the front factor is read straight from our data, accounting for the known entrainment contribution. We now turn our attention to the as yet undetermined time scales τ L and τ B,Loop . τ L is a representative time scale for the loop trajectories, which comes from an average over impact parameters. We do not know its dependence on a p , but we do know it is independent of temperature, as confirmed by Fig. 5(b) where the data at different temperature change sign at the same time. Appendix C describes how to use this fact to isolate τ B,Loop in measurements of initial decay rates knowing only τ rot . This approach finds proportional to Pe rel as expected, and associated with diffusion over a distance close to A s .
Before progressing, we note that this decorrelation appears as a simple exponential because Brownian motion outside the entrained volume lacks the asymmetry that provided the square root in the analogous factor in the entrainment term.
Finally, we can feed τ B,Loop back into our fitted initial decay rates and solve for τ L . This leads to Fig. 6(b), where, similarly to C (p) vv (0), we find it to be constant below R 1, while it rises with an apparent power law ∼ a 1/2 p above.
While the small R behaviour has the usual explanation that the passive particles are behaving as infinitesimal tracers, the apparent power law is harder to understand as it disagrees with the scaling argument in Appendix B, which predicts τ L ∼ a 1 p . This discrepancy could come from the periodicity of our system, which was not accounted for in our calculations, or could simply be a sign our data does not extend to large enough a p to see the expected behaviour. For the purposes of this work it is sufficient to write down the empirical form as With this, our ansatz form for C

VI. ACTIVE DIFFUSION
We now move from the velocity autocorrelation to the active diffusivity, obtained by integrating C (p) vv as per Eq. (4). This we show for both numerical integration of the simulation data and analytical integration of our ansatz C (p) vv using the approximate expressions for parameters in the previous section. These are shown in Fig. 7.
Beginning with our data plotted against a p /a s in Fig. 7(a), we observe two main features: non-monotonicity with a turning point just below a p /a s = 1; and a temperature-dependence on the small a p side. The non-monotonicity requires different physics to be dominating at different regimes. Using the understanding from previous sections, the presence and absence of entrainment at small and large a p respectively accounts for this behaviour, and is supported by the turning point being close to R = 1.
The quantitative behaviour in the two regimes can also be understood using our analysis of C (p) vv . The decay on the large a p side comes primarily from the decay of C (p) vv (0), which is not fully compensated for by the increase in τ L , at least not over the range of our data.
Note that τ L ≪ τ r&t ≪ τ rot , τ B,Loop in the T = 0.1T 0 data here, meaning the attenuation of  Work by Pushkin and Yeomans 40 has argued that the contribution from far field loops truncated by run and tumble events leads to a constant value of D A independent of the run length. We expect our additional decorrelation mechanisms to fall under the same formalism, and hence would expect a constant term in D A that might be seen if we extended our range of a p /a s . However, feeding our parameters into their calculation would put the value of this constant at D A /a s v s φ s = 20.25, which is clearly missing or greatly reduced in our data. We believe our periodic boundaries are the cause of its absence since it is an effect dominated by flow fields at impact parameters of order the run length. In our case, the run length equals the side length of our box, so there will be significant interference from the swimmer's periodic images.
Our final comment on the large a p regime is that the collapse of the curves at different temperatures here is misleading. As Fig. 5(b) showed, there is a significant difference in C (p) vv here and the dynamics truly are affected by the temperature. We believe our 10T 0 data happened to have an equal loss of the negative tail and initial positive decay, but this is not generally true and we expect an intermediate temperature, e.g. 5T 0 , would have a higher D A here because it's negative tail will have been affected the most by the decorrelations. By the same reasoning, a temperature larger than 10T 0 would have smaller D A because there is very little of the negative tail left to lose, leading to a greater loss from the positive part.
The variation at small a p is driven primarily by diffusive processes, so our C (p) vv analysis predicts this behaviour to be a function of τ rot ∼ Pe s ∼ T 0 /T and τ B,Ent , τ B,Loop ∼ Pe rel , instead of a function of a p /a s . In practice we find τ rot is large enough that it has a negligible influence, leading to our data falling onto a master curve when plotted against Pe rel , as shown in Fig. 7(b). The apparent power law we observe across our data is close to Pe 0.2 rel , not Pe 1 rel or Pe 1/2 rel as we might have expected from the form of the two Brownian decorrelation factors.
Our ansatz function provides and explanation for this, although we must first specify a value of c 2 . The first value we use in Fig. 7(b) is c 2 = 0.22, which was chosen by a least squares fit of the ansatz to the data, setting all other parameters as described in Section V A.
Here we see good agreement with the data at all temperatures, including the appearance of a shallow apparent power law. Extending the plot down to smaller Pe rel finds the expected Pe 1 rel behaviour does appear eventually. Importantly, the ansatz plots have undulations, which are made extreme when using the largest allowed value of c 2 = 0.5. This results from the two terms in the ansatz, with the loop term providing the peak at small Pe rel and the entrainment providing the second rise. These undulations are more subtle in our simulation data, but they are still visible in the curvatures of the 10T 0 and 1T 0 data. Hence we attribute the small power law to a transitional regime between loop and entrainment dominance.
It is useful here to compare to Kasyap, Koch and Wu's calculation of D A ∼ Pe 1/2 for small Pe in a slender-body swimmer model 32 . In contrast, both terms in our anstaz lead to D A ∼ Pe 1 in this limit. We suspect the origin of this discrepancy might lie in their result assuming the distance swum in a single run is much smaller than the typical displacement by Brownian motion in the same time, whereas the reverse was true for all systems we used to construct our ansatz.
Finally, we note that the properties of our swimmers, especially the lack of steric interactions, will limit the applicability of our understanding to experimental systems. Our nono-monotonicity is nevertheless in qualitative agreement with the experiments of Patteson et al. 33 This encourages us to propose that the cause is to be found in the transition between a regime dominated by entrainment events for small passive particles, and one for larger particles where far field loops are most important. We believe that this prediction could be testable with the experimental trajectories already available from the experiments in Ref. 33.

VII. CONCLUSIONS
We have looked at the effect of both particle size and temperature on the active diffusion of spherical passive particles in 3D periodic systems of microswimmers. For this we used a 'smart' version of the Wavelet Monte Carlo dynamics algorithm to simulate active systems with hydrodynamically correlated rotations and translations, biased by swimming force. This gave us an efficient algorithm that includes correlated thermally-driven Brownian motion that is sensitive to particle size.
Our first results were geared towards validating active, non-thermal behaviour in smart WMCD, for which we simulated the trajectories of single passive particles at varying impact parameters from a passing swimmer at very large Péclet number. These results were consistent with previous work, demonstrating the expected cusped-loop trajectories at large impact parameter, whose net and maximum displacements decayed with the expected power laws.
We then turned our attention to dilute mixtures of swimmers and passive particles with thermal fluctuations present. By using a range of temperatures and passive particle sizes we were able to identify the physics driving active diffusion via the behaviour of the velocity autocorrelation. Analysis of this led to constructing an ansatz function to unify the diverse forms of C (p) vv observed. This function was expressed as the sum of contributions from entrainment and non-entrained loop trajectories, both subject to exponential decorrelation factors coming from swimmer rotations and Brownian motion. More generally, any mechanism causing passive particles to cross swimmer-induced stream lines could be included in this way.
Most parameters in our ansatz fall under one of two categories: decorrelation times that vary with the appropriate Péclet number; and parameters describing the Pe → ∞ limit governed by the comparison between the hydrodynamic response of the passive particle and the geometric size of the swimmer. By itself, the second category of parameters leads to D A having two regimes when plotted against a p , with a decay away from the flat, small-a p regime where particles act as infinitesimal tracers. The decorrelation factors then introduce a gradient to the small-a p regime, leading to non-monotonic behaviour.
Plotting D A against Pe rel reveals a master curve for the small a p regime. The behaviour of this master curve over the range of Péclet numbers studied is made complicated by the entrainment and loop contributions appearing and plateauing at different values, with their sum leading to a weak apparent power law.
Finally, we note that the simplicity of our swimmer model means it is not expected to give quantitatively relevant results for comparison with experiment. Instead, the strength of our results lies in the identification of the role of particle size and Péclet number(s) in the velocity autocorrelation. In the process we highlighted the importance of temperature and near-field effects, both of which are often neglected in theoretical and computational studies of similar systems.

ACKNOWLEDGMENTS
We gratefully acknowledge funding by the EPSRC, grant no. EP/M508184/1, and the Warwick SCRTP for computational resources. We also thank M. Polin for several discussions and his help with preparing the manuscript.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Appendix A: Loop size calculations
We first consider the noise-free transits as graphed in Fig. 2. The displacement of a passive particle ∆(t) due to the passage of a swimmer incident with impact parameter ρ and swimmer velocity v sẑ obeys where v(r) is the flow field established by the swimmer andρ is the radial unit vector in cylindrical polar coordinates.
As constructed, the swimmer approaches from below so that ∆ z − z s decreases from +∞, down to −∞ for an orbit which does not get entrained. It is then convenient to write v(r) = v s g(ρρ, z − z s ) and v s dt = dz s so z s is the vertical rise of the swimmer, which leads For the unentrained trajectories we can follow earlier work 49 in expanding this for the total deflection as a series in g . It is convenient then to parameterise in terms of z = −z s which is to zero'th order the height of the passive above the swimmer, giving Expanding ∆ in implied powers of g then gives ∆ = ∆ (1) + ∆ (2) + O(g 3 ) where and On the LHS ρ and z parameterise the transit in terms of impact parameter and time (as −z/v s ) through it. However on the RHS they are simply cylindrical polar coordinates of the flow around the swimmer.
We now focus on swimmers with azimuthal symmetry, so we write which can be thought of as a divergence free flow field. Moreover, in the far field the swimmer flow is proportional to that of a force dipole −∂(G TT O ·ẑ)/∂z so we infer that in the far field h approaches G TT O ·ẑ, where G TT O is the Oseen tensor. It then follows that for any force free swimmer the first order advective deflections of a point passive particle form a closed loop ending up with ∆ (1) (ρ, −∞) = 0.
At second order we now need forces separated by A s and with dipole strength κ = A s F s , near approaches ρ ≪ A s lead to ∞ −∞ h ρ (ρ, z ′ ) 2 dz ′ ∝ +ρ/A 2 s and hence ∆ (2) z ∝ +1/(ρA 2 s ), this time with a positive sign. The full forms for a swimmer modelled by two point forces separated by A s can be found as follows. The transverse component of the Oseen tensor is given by (ρz/8πη)(ρ 2 + z 2 ) −3/2 and integrating this with respect to z gives −(ρ/8πη)(ρ 2 + z 2 ) −1/2 , which then leads to where ζ = z/ρ and α = A s /2ρ. In the limit of large α the integral over h ρ (ρ, z ′ ) 2 is then dominated by two well separated Lorentzians each of width ρ and height ∝ A −2 s , leading to ∞ −∞ h ρ (ρ, z ′ ) 2 dz ′ ∝ ρ/A 2 s as used in the paragraph above. To get the full result we write where K(x) = π/2 0 dθ(1 − x 2 sin 2 θ) −1/2 is the complete elliptic integral of the first kind.
in our simulations, measuring the initial decay rate of C vv gives the full decay rate Note we included the quadratic term in the fit, with c 2 /τ 2 L as a second fitted variable, although we do not use those values in this work.
The different terms in τ tot are all expected to have different behaviours across our simulations: τ r&t is a constant; with fixed swimmer parameters, τ rot depends only on the temperature; and τ L depends on particle sizes but not the temperature. By looking at the difference of τ tot at different temperatures but the same a p , we can remove the influence of both τ r&t and τ L . We then only need to remove τ rot , which can be done by hand since we know its form, leading us to consider Assuming a power law τ B,Loop (a p , T ) ∝ Pe α rel we would have so we should still see the same power law if we are consistent with which temperature we use in Pe rel . We have 3 combinations of temperature: 10T 0 and 1T 0 ; 10T 0 and 0.1T 0 ; and 1T 0 and 0.1T 0 . For all combinations we use the smaller temperature in Pe rel , giving the results in Fig. 8. Here we see T ∼ Pe 1 rel , so hence so does τ B,Loop , which is consistent with other time scales increasing linearly with a Péclet number. Furthermore, Fig. 8 provides us with the constant of proportionality, which is 9 times larger than the one for T due to the factor of (T i /T j ) 1 − 1. Hence, we have Errors shown are 1 standard deviation, and are smaller than the plot markers for many data points.