Collisions of particles advected in random flows

We consider collisions of particles advected in a fluid. As already pointed out by Smoluchowski [Z. f. physik. Chemie XCII, 129-168, (1917)], macroscopic motion of the fluid can significantly enhance the frequency of collisions between the suspended particles. This effect was invoked by Saffman and Turner [J. Fluid Mech. 1, 16-30, (1956)] to estimate collision rates of small water droplets in turbulent rain clouds, the macroscopic motion being caused by turbulence. Here we show that the Saffman-Turner theory is unsatisfactory because it describes an initial transient only. The reason for this failure is that the local flow in the vicinity of a particle is treated as if it were a steady hyperbolic flow, whereas in reality it must fluctuate. We derive exact expressions for the steady-state collision rate for particles suspended in rapidly fluctuating random flows and compute how this steady state is approached. For incompressible flows, the Saffman-Turner expression is an upper bound.


Introduction
Turbulent aerosols are of interest in a variety of natural and technological systems. Two very important examples are water droplets in turbulent rain clouds [1] and dust grains in turbulent accretion disks around growing stars [2]. In both of these systems the aerosol is unstable because the suspended particles collide (leading to aggregation or possibly fragmentation of the aerosol particles). The collision processes therefore have significant consequences: the formation of rain in one case, and the widely hypothesised mechanism for the formation of planets in the other. Collisions always occur due to molecular diffusion, but (as pointed out by Smoluchowski [3,4]), macroscopic fluid motion can considerably increase the collision rate.
If the suspended particles are sufficiently heavy (so that their inertia becomes relevant), they can move relative to the fluid. In this case the occurrence of 'caustics' will typically increase the collision rate by several orders of magnitude [5,6]. If the aerosol particles are sufficiently light, their 'molecular' diffusion can make a significant contribution to the collision rate, which can be estimated by standard kinetic theory.
In this paper we are concerned with the effect of macroscopic motion of a fluid on small particles which have insignificant inertia, so that they follow the flow (advective motion). The seminal papers in this area were due to Smoluchowski [3], who first considered the effect of shear of the fluid flow on collisions, and Saffman & Turner [4], who gave a formula for the collision rate which has formed the basis for most subsequent work on this problem. Their paper was motivated by a problem in meteorology argued that small-scale turbulence in convecting clouds can accelerate collisions between microscopic water droplets, thus initiating rain formation: in this case the particles are brought into contact by hyperbolic or shearing motions of the turbulent flow.
The formula for the collision rate by Saffman & Turner [4] has been used frequently in the past five decades in cloud physics and in chemical engineering problems. It appears to be widely accepted that their expression is an exact relation for the collision rate in a dilute suspension. In the following we show that the Saffman-Turner estimate describes an initial transient of the problem only. For particles suspended in incompressible flows, the collision rate falls below the initial transient (which thus constitutes an upper bound). For particles advected in a compressible flows, however, homogeneously distributed particles will cluster in a compressible fluid (see for example [7,8]). The clustering may increase the collision rate beyond its initial transient.
The Saffman-Turner approximation treats the flow surrounding a test particle as if it were a steady hyperbolic flow, while in reality the flow fluctuates as a function of time. In section 3.2 below we give an extension of the Saffman-Turner formula which does give the collision rate exactly. Unfortunately the formula contains information about the time-dependence of the flow, and it is impossible to evaluate it in the general case.
Because of the importance of understanding collision rates for aerosol particles, it is desirable to find exactly solvable cases which can be used as a benchmark for numerical studies. The collision rate must depend on a dimensionless parameter describing how quickly the fluid velocity fluctuates, the 'Kubo number' [9]. We are able to obtain precise asymptotic results on the collision rate in the limit where the Kubo number approaches zero. In this case, particle separations undergo a diffusion process. By solving the corresponding Fokker-Planck equation we can determine the collision rate exactly.
The remainder of this paper is organised as follows. In section 2 we introduce the equations of motion and the dimensionless parameters of the problem. Section 3 discusses the Saffman-Turner theory and our extension of it. Section 3.1 describes the expression for the collision rate given in [4] which is the starting point of our discussions. (Some new results on the evaluation of the Saffman-Turner expression are described in the appendix). In section 3.2 we discuss our exact formula for the collision rate, and explain why the Saffman-Turner approximation describes an initial transient only. In sections 4-6 we discuss how exact asymptotic results may be found for the limit of small Kubo number. A Fokker-Planck equation for the probability density of the separation of particles is described in section 4. In section 5 this is used to obtain the steady-state collision rate and section 6 gives the full time dependence of the collision rate. These results are compared to numerical simulations in section 7, which also contains some concluding remarks, discussing scope for further work in this area.
Finally, we remark that a brief summary of some of the results of this paper has already been published [10]. Here we discuss the problem in greater depth and generality, and derive expressions for the time-dependent collision rate which were not discussed in [10].

Equations of motion and dimensionless variables
We consider spherical particles of radius a in a fluid with velocity field u(r, t) which has an apparently random motion, usually as a result of turbulence. We assume that the suspended particles do not modify the surrounding flow. When the inertia of the particles is negligible, they are advected by the flow: It is assumed that direct interactions between the particles can be neglected until they collide. In other words, the particles follow equation (1) until their separation falls below 2a.
We model the complex flow of a turbulent fluid by a random velocity field u(r, t). We consider flows in both two and three spatial dimensions and for convenience we use a Gaussian distributed field when we carry through concrete computations. In most cases we are concerned with incompressible flow, satisfying ∇ · u = 0. Particles floating on the surface of a fluid may experience a partly compressible flow [11,12], as may particles in gases moving with speeds comparable to the speed of sound. For these reasons we also consider partially compressible flows.
It is convenient to construct the random velocity field u(r, t) from scalar stream functions or potentials [9]. In two spatial dimensions we write u = N 2 (∇ ∧ ψn z + β∇φ) (2) where N 2 is a normalisation factor and ψ and φ are independent Gaussian random functions. We shall use angle brackets to denote averaging throughout. The fields φ(r, t) and ψ(r, t) have zero averages, φ = 0, ψ = 0 and they both have same correlation function, C(R, T ): where R = |r − r ′ | and T = |t − t ′ |. This two-point correlation function C(R, T ) is a smooth function decaying (sufficiently rapidly) to zero for large values of R and T . For β = 0, the flow (2) is incompressible. For finite values of β it acquires a compressible component. In the limit of β → ∞ the flow is purely potential. In some cases the physics of a problem dictates that ψ and φ should have different correlation functions; many of our results can be generalised in this way.
In three spatial dimensions we write where A = (A 1 , A 2 , A 3 ) and φ are four independent scalar fields with zero mean and with the same correlation function C(R, T ); and N 3 is a normalisation factor. In the remainder of this paper we choose the normalisation factors to be of the form where u 0 is the standard deviation of the magnitude of the velocity and C ′′ denotes the second derivative of the correlation function (3) with respect to its first argument. Fully-developed turbulent flows have a power-law spectrum in the inertial range, covering a wide band of wavenumbers [13,8]. This feature can be incorporated by giving C(R, T ) a suitable algebraic behaviour over a range of values of R, as explained in [10]. The long-ranged behaviour of the velocity field is not, however, relevant to the advective collision mechanism. It therefore suffices to consider a model with a short-ranged velocity correlation: for the numerical work reported in this paper we used following form of the correlation function where C 0 is a constant. In our numerical simulations we represent the flow field by its Fourier components which are subject to an Ornstein-Uhlenbeck process as suggested by Sigurgeirson & Stuart [14]. Our problem is characterised by the six parameters listed in table 1. From the parameters in table 1, three independent dimensionless combinations can be formed: Parameter Symbol Particle size a Typical velocity fluctuation u 0 Correlation length of the flow η Correlation time of the flow τ Number density of particles n 0 Compressibility β Spatial dimension d Table 1. Parameters of the model.
The first parameter characterises the dimensionless speed of the flow and is called Kubo number, discussed in [9]. Note that Ku ≫ 1 is not possible, because the motion of the fluid places an upper limit on the correlation time. Steady, fully developed turbulence corresponds to Ku ∼ 1. The Kubo number can be small for randomly stirred fluids. It is only in the limit Ku → 0 that we are able to obtain precise and explicit estimates for the collision rate: this case is considered in sections 4-6. The second parameter in (7) is the packing fraction of particles. Throughout we assume that this parameter is small. Similarly, the third parameter in (7) is usually taken to be small. We note that Kalda [15] has considered the collision rate for particles in a nonsmooth velocity field: this could be relevant to the case where a ≫ η.

The Saffman-Turner expression
Throughout this paper we consider the rate of collision of a given particle with any other particle, denoting this by R. Some papers consider the total rate of collision per unit volume. If the spatial density of particles is n 0 , the total rate of collision per unit volume is 1 2 n 0 R (the factor of 1 2 avoids double-counting). We will not be concerned with what happens after particles undergo their first collision: in different physical circumstances they may coalesce, scatter, or fragment, but in this paper we are concerned only with their first contact.
We assume that the particles are spherical (or circular, in two-dimensional calculations) and that they all have the same radius, a. We regard the particles as having collided when their separation reaches 2a, and we neglect effects due to the interaction of the particles and the fluid. In practice the fluid trapped between approaching particles may slightly reduce the collision rate [4], but this effect can be accounted for by replacing the radius by an effective radius.
The problem of calculating the rate of collision therefore reduces to the following problem. We consider a given particle, and transform to a frame where the centre of this particle is at the origin, and the separation of the centre of another particle is denoted by R. Initially, the reference particle is surrounded by a gas of particles with spatial density n 0 . We assume that these are initially randomly distributed, apart from the constraint that none of the particles is in contact with the reference particle. Collisions with the test particle occur when other particles come within a radius 2a of the reference particle. The rate of collisions is therefore the rate at which particles cross a sphere of radius 2a centred at the origin in the relative coordinate system. This is obtained by integrating the inward radial velocity over the sphere, and multiplying by the density n 0 . This approach gives an expression for the collision rate which we term R 0 : where v r (R, Ω, t) is the radial velocity at spherical coordinate Ω, radius R and time t. The function Θ(x) is a Heaviside step function, which is used to select regions of the surface where the flow is into the sphere of radius 2a. This is the fundamental expression for the collision rate given by Saffman & Turner [4]. In section 3.2 we discuss why this expression is not exact, and give the precise formula. The remainder of this section considers how this expression is evaluated; some new results are presented in the appendix.
The evaluation of (8) is greatly simplified in the case where the particles are small, in the sense that a/η ≪ 1. In this case the relative velocity is accurately approximated using the velocity gradient of the random velocity field u. This approximation was also considered by Saffman & Turner [4]. To lowest order in R, the relative velocitẏ R = u(r + R, t) − u(r, t) is approximated as A(r, t)R where A is the rate of strain matrix of the flow u, with elements A ij = ∂u i /∂r j (and i, j = 1, . . . , d). Two particles with radii a ≪ η moving close to each other will thus experience a relative velocity A 0 R. This is illustrated in figure 1 for a flow which is hyperbolic in the vicinity of the reference particle. The corresponding relative radial speed v r is v r =n T A 0 R wheren is the radial unit vector. In particular, at the distance R = 2a where particles collide, the radial speed is v r = 2an T A 0n . Thus when a/η ≪ 1, equation (8) reduces to This expression was also obtained in [4]. The remainder of this subsection is concerned with the evaluation of (9). Earlier, Smoluchowski [3] had considered the special case of a fluid flowing with a uniform shear in 3 dimensions: He obtained the collision rate: which is in agreement with the result of evaluating (9) for this case. For a general strain-rate matrix A 0 the evaluation of (9) is, however, very difficult. In Appendix A we discuss how this expression is evaluated for a general matrix A 0 in two dimensions, and for a general traceless matrix (representing an incompressible flow) in three dimensions. In a turbulent flow, Saffman & Turner [4] argued that the elements of A 0 change as a function of position and one needs to average over the ensemble of strain matrices A 0 at different positions in order to estimate the collision rate, so that (9) is replaced by: At first sight the requirement to average over A 0 appears to complicate the problem. However, for a rotationally invariant ensemble of random flows, the problem is considerably simplified by taking the average. In an incompressible flow, for each realisation of u(r, t), the currents into and out of the collision region (the disk or sphere of radius 2a) cancel precisely. Therefore (12) can be written as The same result also holds for cases where the flow is compressible, of the form (2) or (4), with β = 0. This is demonstrated by the following argument (here we discuss the two-dimensional case). If we did not include the factor Θ(−n T A 0n ) in (13), we would calculate the sum of the collision rate for the flow u and for the time-reversed flow −u.
The time reversed flow is generated by reversing the signs of the potentials φ and ψ in (2). The probability density for the Gaussian field −φ is the same as for φ (and similarly for ψ). It follows that the collision rate for the time-reversed flow −u is the same as for u. The expression (13) is therefore also valid for the cases of compressible flow which we consider in this paper. Now using rotational symmetry, one finds the very simple expression where A d (2a) is the the area of the sphere of radius 2a in d dimensions (explicitly A 2 (r) = πr 2 and A 3 (r) = 4πr 3 /3). For the Gaussian model flow which we consider, equation (14) gives with

An exact expression for the collision rate
The Saffman-Turner expression for the collision rate, equation (8) or (14), correctly describes the initial collision rate. There are, however, two reasons why the collision rate may approach a significantly different value after an initial transient. The first reason why (8) may fail arises from the fact that the flow field fluctuates in time. Recall that we are concerned with the rate at which pairs of particles collide for the first time. If the flow is time-dependent, the relative position coordinate R(t) may pass through the sphere of radius 2a more than once. This effect may be accounted for by writing the collision rate in the form As before dΩ is the d-dimensional surface element at R = 2a and v r is the radial velocity component (the relative speed) and the Heaviside step function ensures that only particles entering the sphere contribute to the collision rate. The factor n(t) is the density of particles in the neighbourhood of the test particle at time t. The function χ is an indicator function: it is equal to unity if the point reaching the surface element Ω at radius 2a and at time t has not previously passed through the sphere of radius 2a, otherwise it is zero. The effect of including the function χ is illustrated in figure 2.
The collision rate will reduce below the Saffman-Turner estimate for times where χ is no longer unity.. In the general case (17) is difficult to evaluate since the indicator function χ depends on the history of the flow. If the flow is rapidly fluctuating (that is, if the Kubo number of the flow is small), the relative separation of two particles undergoes a diffusion process which makes it possible to exactly evaluate the collision rate for a given ensemble of random flows.
A second effect may modify the collision rate from the Saffman-Turner estimate is that if the flow field is compressible, particles may cluster together, and the particle density in the vicinity of a test particle may be higher than the average density. This effect is expected to cause the collision rate to increase, after an initial transient during which the clustering becomes established. Again, this effect is very hard to quantify in the general case, but we can obtain precise asymptotic results in the limit of small Kubo number, where a diffusion approximation becomes applicable. Because we consider particles which are advected by the flow, there is no clustering effect when the flow is incompressible.
We note that for incompressible flow the only correction to the Saffman-Turner estimate is the occurrence of the factor χ in (17). This implies that for advected particles in an incompressible flow, the Saffman-Turner estimate of the collision rate is an upper bound.

Diffusion approximation for small Kubo number
In the limit of a rapidly changing flow, that is when Ku ≪ 1, particle separations R = r−r ′ undergo a diffusion process (see [8] for a review). By solving the corresponding Fokker-Planck equation with the appropriate boundary conditions, we can determine the collision rate exactly in this limit.

Fokker-Planck equation
Consider two particles, one at r and the other in its vicinity, at r ′ . The equation of motion (1) implies that their separation R = r ′ − r obeyṡ Integrating (18) over a short time interval δt we obtain The moments of the components δR i of δR are: assuming that δt ≫ τ . We thus obtain the following Fokker-Planck equation for the density ρ(R, t) of particle separations [7,8,10]: Here D(R) is a diffusion matrix with diffusion coefficients In terms of the correlation function C(R, T ) introduced in section 2 we have In order to determine the collision rate it is convenient to write the Fokker-Planck equation in the form of a continuity equation ∂ρ ∂t where the components of the probability current can be identified as j = −D(R)∇ρ. The collision rate R(t) between particles of radius a is given by the rate at which the particle separation R decreases below 2a [10]. This rate equals the radial probability current of particle separations evaluated at R = 2a, i.e.
wheren is the radial unit vector. We evaluate (25) in two steps. First (in section 5) we determine the steady-state collision rate as t → ∞. Second, we determine the full time-dependence of the collision rate (section 6), describing how the initial transient discussed in section 3 approaches the steady state.

Transformation to spherical coordinates
Because of the angular symmetry of the statistics of the fluid velocity, it is convenient to transform (21) to spherical coordinates. The probability density ρ depends on the separation R only and obeys Here the radial probability current j r (R, t) is given by with

Steady-state solution
Consider now the steady-state solution of (21), obeying ∂ρ/∂t = 0. The steady-state collision rate, denoted by R ∞ , is determined by (25) in terms of the steady-state current. The total current −A d (R)j r (R) entering a sphere of radius R and area A d (R) must be a constant. This gives j r (R) = −R ∞ /A d (R) where R ∞ is a constant, equal to the collision rate which we wish to determine. The equation of motion (1) is only physically meaningful when a ≪ η, so that the fluid velocity is approximately the same throughout all of the region occupied by the particle. We will, however, solve the diffusion equation for general a. One reason for treating the general case is that it is hard to verify the analytical expressions for the limiting case a/η ≪ 1 in numerical investigations.
To determine this current, it is necessary to consider the appropriate boundary conditions. First, if the particle distribution initially is uniform with density n 0 , we expect it to remain so at large separations. Thus we use the boundary condition The boundary condition (29) is implemented as follows. For a specified current density j r (R, t) we can solve (27) by finding an integrating factor h: with integration constant A and The boundary condition (29) determines the integration constant in (30) to be A 0 . Second, when particles comes closer than the distance of 2a they collide and must be removed. This is taken care of by the following boundary condition at 2a: Inserting this boundary condition into (30) yields Using j r (R ′ ) = −R ∞ /A d (R ′ ) and solving for the constant R ∞ we find the following expression for the collision rate: where Γ(z) is the Gamma function. Equation (34) is the main result of this section. We continue by discussing a number of limiting cases.

Collision rate for incompressible flows
For an incompressible flow (β = 0), the integrating factor is just unity (incompressibility ∇ T u = 0 implies ∇ T D = 0 which in turn implies that the term proportional to ρ in (27) vanishes). In this case, equation (34) simplifies to

Collision rate for β = 1
When β = 1, the strength of the solenoidal and potential parts of the field are equal. This case may be relevant to the dynamics of particles floating on the surface of a turbulent fluid [11]. The collision rate can be evaluated exactly in this case, by rewriting (34) as Setting β = 1 we find for d ≤ 2. When d > 2 we obtain What makes it possible to find exact solutions in this case is the fact that when β = 1, the functions f (R) and g(R) in (28) are related as g ′ (R) = f (R). Note that this relation is always true in one spatial dimension, independent of the value of β. Note also that when g ′ (R) = f (R), equation (27) can be solved for ρ by integration (starting at e.g. 2a, with ρ(2a, t) = 0) Thus, a non-vanishing steady-state collision rate is obtained only if g ′ (R) = f (R) or when d > 2.

Collision rate for small particles
The collision rate must vanish as the particle size approaches zero. This means that the integral in the denominator in (34) must diverge for small a. Thus if a ≪ η is small enough, the major contribution in the integral in (34) comes from small values of R and the relative error will be small if we replace the integrand by its small R expansion. It is convenient to change variables according toR = R/η. We expand and obtain approximately Here the parameter D is defined as It corresponds to the diffusion constant D 1 /(mγ) 2 in equation (39) of [9] and in the incompressible case it corresponds to the diffusion constant D defined in equation (17) of [10]: here Γ is an alternative parametrisation of the compressibility, defined as follows The parameter Γ was also used in [9] (and earlier work cited therein). It can take values between Γ min = 1/3 and Γ max = (d + 1)/(d − 1). Here Γ max corresponds to a completely incompressible flow with β = 0 and Γ min corresponds to a purely potential (β → ∞) compressible flow. In particular, when Γ = 1 the field strengths of the compressible and incompressible components of the flow are equal, β = 1. Using equations (41), we get the collision rate for small particles This expression was previously derived in [10].
For an incompressible flow, with Γ = (d + 1)/(d − 1), the particle density is uniform and the collision rate is proportional to the packing fraction n 0 a d as expected. By contrast, if the flow is compressible, the particles will cluster on a fractal set [16]. This is expected to enhance the collision rate because the particle density is large within the clusters. The clustering effect can be characterised by the correlation dimension of the particles, D 2 . It is defined by the scaling law P (ǫ) ∝ ǫ D 2 , where P (ǫ) is the probability that the particle separation is smaller than ǫ [17]. We have and thus D 2 = (d − 1)Γ − 1. This expression was derived in [7]. Equation (44) shows that in compressible flows, the collision rate depends upon the correlation dimension D 2 = (d − 1)Γ − 1 < d. Whenā ≪ 1 the collision rate in a compressible flow is therefore much larger than the corresponding rate in an incompressible flow.
Note finally that the steady-state collision rate (44) tends to zero when Γ → For values of β larger than β c , no steady-state current satisfying our boundary conditions exists. In one spatial dimension the parameter β c vanishes: no non-trivial steady state exists because all particle trajectories eventually coalesce (this effect was termed 'path-coalescence' Wilkinson & Mehlig [18]).

Collision rate for the Gaussian correlation function
When the correlation function C(R, T ) is of the Gaussian form (6), we obtain from (34) where as before R = ηR. The major contribution to the integral comes from small values ofR (except for when d = 2 and β = 1 or d = 1, when the integral diverges) and for small values ofā, the integrand can be expanded in powers ofR. Expanding the exponential function in the integrand yields where the function h was defined in (31). The diffusion constant D is given by (42). With (6) we find For low dimensions, is a fairly good approximation (a maximum of 3 percent error when d = 2). Substituting (49) into (47) gives once more (44).

Time-dependent solution
So far we have derived the probability densities and collision rates for the steady state.
We now wish to derive expressions of these quantities as functions of time. To this end we need to solve the time dependent Fokker-Planck equation (26) with the boundary conditions ρ(2a, t) = 0 , and ρ(∞, t) = n 0 .
As before we assume a uniform initial scatter of particles The collision rate is given by equation (25) where we have used that the area of a d-dimensional sphere of radius r is A d (r) = 2π d/2 r d−1 /Γ(d/2), and the gamma function is denoted by Γ(z). In order to render the boundary conditions homogeneous, we split ρ(R, t) = ρ 1 (R, t) + ρ 2 (R) and impose the conditions ρ 1 (2a, t) = ρ 1 (∞, t) = 0, and that the differential equation homogeneous in ρ is also homogeneous in ρ 1 . Thus ρ 2 (R) is uniquely given by the steady-state solutions found in the previous section. Now consider the remaining equation for ρ 1 (R, t) which is identical to (26) with homogeneous boundary conditions and initial condition Separation of variables ρ 1 (R, t) = F (R)G(t) in (26) gives (using the radial current (27)) Since the left-hand side depends on t only, and the right-hand side depends on R only, both sides must be equal to a constant, B say. Choosing B = −Dµ 2 , where µ is a positive dimensionless constant, we obtain Note that the solution corresponding to µ = 0 has already been taken care of in ρ 2 . Therefore we can restrict ourselves to considering µ > 0 in the following.
To solve the radial part of (54), we consider the limit of R ≪ η, where f (R) and g(R) are simple power laws (this follows from (41)). Using the dimensionless variablē R = R/η gives the following approximate equation for small values ofR This is an Euler equation. Its solution is obtained by the variable substitutionR =āe u , where the factorā = 2a/η is included in for later convenience and 0 ≤ u < ∞. We find where and D ± and φ ± are integration constants.
We insert this expression into the ansatz for ρ 1 . Upon changing order of the integrations (that is we perform the s-integral first), we find To calculate the collision rate we need to know the current j r at radius 2a. To this end we just require the derivative of ρ evaluated at R = 2a (since ρ itself is constrained to vanish there). We find Equation (52) now allows us to calculate the time-dependent collision rate. We split j r [see equation (27)] into two parts depending on ρ 1 and ρ 2 respectively. The second part gives the steady-state collision rate R ∞ , see equation (34). We find (expanding in powers ofR) This is our main result for the time-dependent collision rate, valid for correlations with non-vanishing fourth order derivative at R = 0 and for small particles. An approximation to this result can be obtained by approximating the steady-state probability density by expanding in powers ofR. Since theũ-integral extends to infinity and since the Gaussian contribution in the integrand becomes small for large values of Dt, we must have that ρ(ũ, ∞) goes to n 0 asũ goes to ∞, which is the case for the original ρ(ũ, ∞).
To accomplish this, we match the small-R expression of ρ(ũ, ∞) atR = 1 to the large-R expression. By using (41) and (44) and matching atR = 1 we find where the last term is not cut off atR = 1 as it vanishes sufficiently quickly. The time-dependent collision rate can now be evaluated from the sum of three integrals of the type Putting everything together we find valid forā ≪ 1. For times less than [d(d + 2µ 0 )D] −1 the above expression simplifies to For intermediate times, while for large times, t ≫ ln(ā)/(2µ 0 D), where Consider finally evaluating (71) in the incompressible limit, we have µ 0 = −d/2 and the terms involving the integral cutoffs ln(ā) cancel (I 1 = −I 2 in (70)). Thus in the case β = 0 the collision rate simplifies to Equations (71) and (75) are compared to results of numerical simulations in the following section.

Numerical illustration and concluding remarks
We performed simulations of the collision rate as a function of time at small Kubo numbers, for both incompressible and compressible flows. These illustrate ( figure 3) the very complex behaviour of the model. For example, in the case of compressible flows the collision rate at first decreases below the value determined by the Saffman-Turner approximation due to the effect illustrated in figure 2, before rising again as particle clustering becomes apparent. PSfrag replacements Figure 3. Collision rate for particles advected in a two-dimensional flow of the form (2) with (6) and C 0 = u 2 0 η 2 /2. Particles are initially randomly distributed, initially overlapping particle pairs are removed, as are particles which have collided. Parameters: n 0 = 1000, u 0 = 1, η = 0.1 and a = 0.001 for all graphs and for a β = 0 and τ = 0.004, b β = 0 and τ = 0.0008, c β = 1/ √ 13 and τ = 0.004, d β = 1/ √ 13 and τ = 0.0008, e β = 1/ √ 5 and τ = 0.0008, f β = 3/7 and τ = 0.0008. The collision rate is approximated by the cumulative sum of all collisions up to t, divided by t. The light blue areas are the intervals R ± 2σ/ √ N , where σ is the standard deviation of the rate, and N is the number of realisations of the flow. The latter decreases over time (it is typically N ∼ 10 4 for the first two decades in t, and N ∼ 10 for the last decade). The Saffman-Turner estimates (15,16) are shown as red dashed lines. Our own result (71) is shown as red solid lines. To correspond to the simulated collision rate, the plotted theoretical collision rates has been integrated to and then divided by t. Due to the assumptions in the Fokker-Planck theory, the long time rate (71) is not valid for 7.1. Simulations Figure 3 summarises our results for the collision rate of particles advected in flows with small Kubo numbers. Shown are numerical simulations for a two-dimensional random flow of the form (2) with correlation function (6).
Panels a and b show the collision rate in an incompressible flow (β = 0). As expected (see section 3.2), the collision rates drops below the initial transient given by the Saffman-Turner approximation, (15,16) with β = 0. For the particular choice (6) with C 0 = u 2 0 η 2 /2, equations (15,16) become in two spatial dimensions Also shown is our own theory valid for small Kubo numbers (equation (71) in section 6). The agreement between the theory and the simulations is good in all cases, but slightly better for the smaller Kubo number (panel b). As discussed in section 3.2, the initial transient constitutes an upper bound to the collision rate.
Panels c to f in figure 3 show the collision rate in compressible flow. Now, the Saffman & Turner theory is no longer an upper bound (see for example figure 3f), because initially homogeneously distributed particles in an incompressible flow cluster together. The corresponding density fluctuations increase the collision rate, as our exact result shows (equation (71) in section 6). Again we observe good agreement between the simulations and our analytical result.

Scope for further investigations
In this paper we have concentrated upon the solvable case of advective collisions in flows with small Kubo number, which provides considerable physical insight. We conclude by commenting on the relation between these results and collisions in a turbulent flow field which satisfies the Navier-Stokes equation. The standard approach, based upon the Saffman-Turner formula, predicts a collision rate R ∼ n 0 a d τ −1 , where τ is the Kolmogorov timescale [13] of the turbulent flow. The calculation based upon the diffusion equation gives a collision rate R ∼ Ku n 0 a d τ −1 . Given that Ku = O(1) for a turbulent flow, we see that the Saffman-Turner and diffusive expressions are of the same order.
These observations are consistent with the hypthesis that the collision rate for small particles in a turbulent flow is where E is the rate of dissipation per unit mass, ν is the kinematic viscosity, and K d is a universal constant (depending only on the dimension). It would be a valuable addition to the literature on aerosols and suspended particles to determine the value of K 3 from numerical simulations using a Navier-Stokes flow.

Appendix A. Evaluation of collision rate for constant shear
In this appendix we show how to evaluate the expression (9) for a general matrix A with elements A ij (we drop the subscript 0). It is convenient to decompose A into a symmetric part A () and an antisymmetric part The collision rate is independent of the antisymmetric part A [] (because rotations do not contribute to the collision rate). The symmetric part A () can be diagonalised by an orthogonal transformation, A () = OSO T , where S is diagonal with the eigenvalues σ 1 ≤ σ 2 ≤ ... ≤ σ n of A () . In evaluating (9) we may writen T An =n T OSO Tn =n ′T Sn ′ , wheren ′ = O Tn . Sincen ′ is just a rotation ofn which is integrated over all directions, we can replace (9) by Appendix A.1. Two spatial dimensions We now show how to perform the integral in (A.2) in two spatial dimensions. Note that if the flow determined by A is not area preserving, the particle density will change as a function of time. In this case n 0 in (9) must be replaced by At short times we approximate n(t) ≈ n 0 and find where σ + = (σ 1 + σ 2 )/2 and σ − = (σ 2 − σ 1 )/2 > 0. This is the final result in two spatial dimensions, expressed in terms of the eigenvalues σ 1 ≤ σ 2 of the symmetric part A () of the strain matrix.