New insights from comparing statistical theories for inertial particles in turbulence: II. Relative velocities

Part I of this two-part series compared two theories for the radial distribution function (RDF), a statistical measure of the clustering of inertial particles in isotropic turbulence. In Part II, we will contrast three theoretical models for the relative velocities of inertial particles in isotropic turbulence, one by Zaichik et al (2009 New. J. Phys. 11 103018), the second by Pan et al (2010 J. Fluid Mech. 661 73) and the third by Gustavsson et al (2011 Phys. Rev. E. 84 045304). We find that in general they describe the relative velocities in qualitatively similar ways, capturing the influence of the non-local dynamics on the formation of caustics and non-smooth scaling behavior in the dissipation range. We then compare the theoretical predictions with direct numerical simulation data and find that although they capture the qualitative behavior of the data consistently, they differ quantitatively, and we discuss the possible sources of error in each of the theories. Finally, we consider how the Zaichik et al theory predicts that the formation of caustics modifies the form of the RDF, and show that the theory describes the behavior of the RDF for particles whose response time scales with the inertial range time scales of the turbulence.


Introduction
This is the second part of a two-part paper investigating theoretical models for describing the relative motion of inertial particles in isotropic turbulent flows. A central theme and motivation of this work is to improve our understanding of how turbulence impacts the collision rate of particles suspended in the flow. This has important implications for the growth rate of droplets in clouds [1][2][3][4], among many other applications (a more detailed background and motivation has been presented in Part I [5] and will not be repeated here). The collision kernel for a monodisperse population of particles in isotropic turbulence can be written as follows [6,7] where σ is the particle diameter, σ g ( ) is the radial distribution function (RDF) evaluated at contact, ∥ w is the relative velocity of a particle pair along the line of centers, and σ ∥ p w ( ) is the probability density function (PDF) for the relative velocity at contact. Part I [5] focused on the spatial clustering of the inertial particles, characterized by the RDF. The definition of the RDF can be found in Part I. In this second part, we focus on the relative velocity statistics, which are the other major input to the collision kernel.
For inertial particles in isotropic turbulence, the relative velocities between particles can be characterized by the N th order relative velocity structure function tensor, , where τ p is the particle momentum reponse time and τ η is the fluid Kolmogorov time scale). When ≪ St 1 the effect of inertia is weak and S N p scales with r in the same way as S N f . However, recent investigations (e.g. [9][10][11]) have shown that for finite St, S p 2 may deviate significantly from the smooth scaling behavior of S f 2 in the dissipation range (i.e. η ≪ r ), and this has been explained in terms of the formation of 'caustics' (also referred to as the 'sling effect' [12,13] and 'random uncorrelated motion' [14,15], although it should be noted that random uncorrelated motion is really a statistical manifestation of caustics, which are instantaneous events). Caustics, a concept put forth by Wilkinson and Mehlig [16], may be explained as follows. Consider figure 1, which depicts a one-dimensional (for simplicity) phase-space manifold with instantaneous particle positions and velocities (we assume a continuous manifold representing an infinite sample of point particles, consistent with most statistical theories). In the region of space between the caustic lines (the locations of the folds) we may say two things about the particle velocities: (a) particles separated by a small but finite distance will exhibit much larger relative velocities than the underlying fluid relative velocity; and (b) particles with zero separation will exhibit finite relative velocities. The first implies ≫ S S N p N f for small but finite r, and we argue that the second implies for > N 0, although we shall see that this is a point of contention among theories that attempt to describe S N p . Caustics form because unlike fluid particles in incompressible flows whose velocity is uniquely defined by the local fluid velocity field, inertial particles posses a non-local contribution to their dynamics imparted by their finite inertia that allows different particles to arrive to the same position at the same time, but with different velocities reflecting their different path histories.
In this paper, we consider theoretical models for predicting S N p . Two of the theories only describe the N = 2 case and so we focus on S p 2 . As shown in Part I of this paper, S p 2 is especially important because it plays a crucial role in the mechanism giving rise to the spatial clustering of particles.
The paper is organized as follows. In section 2 we review the theories, considering how they account for the formation of caustics and the non-local contribution to the relative velocity dynamics. In section 3, we compare the theories with direct numerical simulation (DNS) data in order to assess their quantitative capabilities. Finally in section 4, we consider how the behavior of the relative velocities influence the behavior of the RDF in both the dissipation and inertial ranges of the turbulence. Conclusions are given in section 5.

Theoretical models for S p 2
The theories and numerical results discussed in this paper are all based on the consideration of the relative motion between two identical inertial particles moving under the Stokes drag force only that do not interact with each other through physical collisions or hydrodynamic interactions and that do not affect the turbulence (i.e. 'one-way coupling'). A monodisperse population of particles is suspended in a turbulent flow that is incompressible and statistically stationary, homogeneous and isotropic. In the Lagrangian frame of reference of a 'primary' particle (fixed at the origin of the Lagrangian reference frame, by definition), the equation for the relative motion of a second 'satellite' particle under Stokes drag forcing is [17] is the relative velocity vector and Δu p is the difference between the fluid velocity evaluated at the positions of the primary and satellite particles. Here and throughout, the superscript 'p' on r and w denotes the time-dependent, Lagrangian variable defined along particle trajectories, and whenever such superscripts are absent (i.e. r and w) the variables denote the 'fixed' reference coordinate system (i.e. either the configuration space r or the phase-space r w [ , ]).
which describes the PDF for particle pair relative position t r ( ) p and relative velocity t w ( ) p in the phase-space r w , . This transport equation for p is closed by use of the Furutsu-Novikov theorem (see [18]) and by approximating Δ t u r ( , ), the fluid relative velocity field, as a spatiotemporally correlated Gaussian field. A more detailed exposition of the ZT can be found in Part I of this paper.
Although the ZT begins with the construction of a PDF equation for the full phase-space dynamics of the particles, in [19,20] it is not solutions to the full PDF equation that are sought, but only the first few moments of the PDF. There are two reasons for this. First, as noted in [19,20], the Gaussian approximation for Δu implies that higher-order moments of p t r w ( , , ) are likely to be in error, especially in the inertial and dissipation ranges, which are known to be characterized by intermittency and PDFs that are highly non-Gaussian [21]. As we shall see, even S p 2 seems to be affected by this approximation. Second, the PDF equation in the ZT is difficult, if not impossible to solve (except in trivial cases), due to the inherent non-Markovian properties that arise from capturing the finite temporal correlations of Δ t u r ( , ), which manifests in the PDF equation as negative eigenvalues in the diffusion tensor [22], giving the equation anti-diffusive properties that make it inherently unstable numerically. For these reasons the PDF equation constructed in the ZT is not solved, but is rather used as a master equation from which the moment equations are constructed.
The transport equations for the first two moments of the PDF for a statistically stationary, isotropic system are is the configuration space PDF (proportional to the RDF, see Part I), are the second and third order particle velocity structure functions. The tensors λ and μ are dispersion tensors whose closed form expressions are given by ⎧  [19,20]), τ r is a Lagrangian time scale given by and τ I is the Lagrangian integral time scale for the flow. Note that in [19,20] λ and μ are made continuous functions of r by specifying curve fits for S f 2 and τ r that connect their various ranges in a continuous way.
The first term on the rhs of (5) represents the contribution to S p 2 coming from particle interactions with the local turbulence and it satisfies The second term on the rhs of (5) involves S p 3 . For ≪ r L I there are two contributions to S p 3 : a local contribution arising because of the skewness of the field Δu [23] and a non-local contribution arising from path history contributions to the dynamics of ( ) t w p . The origin of the non-local contribution may be explained as follows. Particle pairs separated by r that originated from larger separations have experienced (on average) more energetic turbulence than pairs that have come from smaller separations in their path history. Because inertial particles have memory, this asymmetry in the energy along their respective paths will give rise to ; that is, a net inward flux. Notice that this flux exists even if the PDFs of Δu are symmetric. With increasing St, it is this path-history-dependent non-local term that gives rise to caustics.
In [19,20], (5) is closed by introducing a closure approximation for S p 3 . The transport equation for S p 3 (obtained from the equation for p t r w ( , , )) contains the fourth-order particle velocity structure function, S p 4 , which is closed by assuming that it behaves like the fourth order moment of a Gaussian distribution ('quasi-Gaussian' or 'quasi-normal' approximation [24], hereafter denoted by QGA) allowing S p 4 to be expressed in terms of products of S p 2 . With this approximation, the following closed expression for S p 3 is obtained for isotropic turbulence ⎡ As noted in [25], the implications and physical meaning of this QGA were not considered in [19,20]. We now discuss some of their implications. As noted earlier, the PDFs of Δu are strongly non-Gaussian in the dissipation range of the turbulence, and at small St this would also imply strongly non-Gaussian distributions for ( ) t w p . Since the ZT approximates Δu as Gaussian, the QGA for S p 4 introduces no new assumptions for small St. Furthermore, it may appear that since particles become less responsive to the underlying fluid turbulence with increasing St, the PDF for ( ) t w p in the dissipation regime should become increasingly Gaussian with increasing St. However, this is not the case. Indeed, the PDFs for in the dissipation range for larger St may be extremely intermittent [26]. Even though a higher Stokes particle is less capable of responding to the local intermittent structures than a lower Stokes particle, its dynamics are also more non-local. The path history effects discussed earlier give rise to significant skewness and heavy tails in the PDF for in the dissipation range [27]. The QGA for S p 4 in the ZT underestimates the effect of these heavy tails; consequently, while the ZT captures the effects of caustics through the second term on the rhs of (5), the magnitude of those caustic events is substantially underestimated because of the closure invoked for S p 3 . For isotropic turbulence equations (4) and (5) can be cast into spherical coordinates and reduced to a dependence solely upon the separation magnitude r, and S p 2 is given in terms of the longitudinal and transverse components of S p 2 , i.e. ∥ S p 2 and ⊥ S p 2 . The isotropic forms of equations (4) and (5) are given in [20]. Figure 2 shows solutions to these equations for ∥ S p 2 . We note that in [20] the boundary conditions employed correspond to the case of perfectly elastic particle collisions at the collision radius r min . The ZT can in principle describe coalescence and inelastic collisions, however one would need to capture this through the specification of appropriate boundary conditions, which is not trivial. In contrast, such boundary conditions are handled relatively easily when directly solving the underlying PDF equation rather than its moment equations (see [28]). For each of the ZT results shown in this paper, = r 0 min corresponding to the 'ghost' or 'point' particle approximation.
S f 2 , so that only non-local effects (arising from the second term on the rhs of (5)) can cause > ∥ ∥ S S p f 2 2 . Secondly, for St = 1.5 there are two clearly distinct regions for ∥ S p 2 , a power-law region for η > r 0.3 followed by a plateau region at smaller separations resulting in a non-zero intercept for ∥ S p 2 in the limit → r 0. The physical explanation for the transition is as follows. For separations above the critical separation, η r c , the fluid velocity differences Δu are energetic enough to influence the motion of the particles, whereas below this critical separation the fluid velocity differences are no longer sufficiently energetic to appreciably influence the particle motion, causing the particle trajectories below this critical separation to be nearly ballistic. Naturally, the transition occurs at increasing separations with increasing Stokes number and can occur within the inertial range or beyond at sufficiently high values. Figure 2 (b) shows a linear plot of the second-order structure function as a function of the separation distance. At larger separations, < ∥ ∥ S S p f 2 2 , implying the particle relative velocity dynamics are growing increasingly local with increasing separation. This is because Δu becomes statistically independent of separation at the integral scales of the turbulence. Moreover, at large separations, the particles become independent of each other, causing the structure function to approach twice the kinetic energy of the particles, a decreasing function of the particle Stokes number due to the filtering effect of particle inertia [29,30]. Note that all of these trends qualitatively match those found in the DNS study of Salazar and Collins [11].
Concerning the effect of the turbulent Reynolds number ( λ Re , the Taylor micro scale Reynolds number), it is important to recall that the ZT makes a Gaussian approximation for the field Δ ( ) t u r, , hence it does not capture the effects of turbulence intermittency that is known to increase with increasing λ Re . However, the ZT captures the 'classical' λ Re effect, namely that the size of the inertial range, and τ τ η I , increases with increasing λ Re , which the ZT accounts for through its expressions for S f 2 and τ r (see [20] for more details). In particular, the ZT correctly predicts the increase in the relative velocities in the dissipation range at higher Stokes numbers with increasing λ Re , as particles retain more of a memory of their interaction with the larger inertial range scales. This behavior is consistent with that found in DNS [31]. , , but instead directly constructs an expression for S p 2 by using the general integral solution for ( ) t w p , averaging, and then making various approximations to derive a closed integral expression for S p 2 . They developed a model for S p 2 for the general case of a bidisperse The tensor , defined below, accounts for inertial differences of the particle pair ij ij where Ω τ τ = p I 1 1 , , τ τ = z T I and τ T is the fluid Taylor time scale. As we are concerned with monodisperse particles in this study, the tensor  vanishes and will not be considered any further. The tensor  is given by the integral expression where Ψ is the autocorrelation for the fluid velocity difference Δu evaluated along particle pair trajectories and depends upon τ TR and τ R , the Taylor and Lagrangian time scales for Δu evaluated at separation distance R. The PT specifies S f 2 , τ TR and τ R as continuous functions of r (the functions they employ for S f 2 and τ R are in fact the same as those for S f 2 and τ r used in the ZT) whose expressions may be found in [25]. The separation R is a function of time and is defined as 1 2 where φ is the rms separation distance between particle pairs which have the separation r at = = ′ ″ t t 0. This variable enters into the PT as a result of their closure approximation for the fluid velocity autocovariance measured along particle pair trajectories arriving at separation r at = = ′ ″ t t 0. In (12), φ is evaluated at times ⩽0 and is therefore the rms separation backward in time, and they highlight the importance of this feature in understanding the behavior of S p 2 .
However, they note that there is at present no theoretical description or data for the backwardin-time separation of inertial particles in isotropic turbulence and as a result they are forced to approximate φ by its forward-in-time equivalent. As we shall see, this approximation is likely a source of error in the PT predictions compared to DNS data, given that we expect backwardand forward-in-time separation of inertial particles to be quite different (they are different for fluid particles, e.g. [32,33]). Based on the findings in [9] for the forward-in-time separation of inertial particles in isotropic turbulence, PT specifies φ by a ballistic phase followed by a Richardson-Obukhov (RO) phase followed by a diffusive phase and t d is defined as the time at which φ = L I . If ⩾ r L I then the pair motion is governed by the ballistic phase followed by the diffusive phase. Note that φ depends upon S p 2 so that the PT equation for S p 2 must be solved iteratively. According to (12), the degree to which the inertial particle velocities differ from those of the fluid is determined by the weighting 2 . However, as the inertia increases the non-local contribution to the particle relative velocities becomes increasingly important, leading to caustics. The PT predicts the following non-zero intercept for S p 2 at zero separation Since particle pairs were separated in the past (i.e. 0for r = 0 except at = = ′ ″ t t 0) then the integrand in (15) may be non-zero even for r = 0.
Whereas the ZT can account for non-elastic collisions between particles at a finite collision radius through the specification of appropriate boundary conditions, the PT implicitly assumes elastic ghost collisions. It is likely that φ is the main parameter in the PT that would have to be modified in a non-trivial manner in order to account for the effect of inelastic collisions.
In figure 3, we plot solutions to (12) for ∥ S p 2 for monodisperse particles. While there are obvious quantitative differences between the predictions of the PT and the ZT in the dissipation regime (cf figure 2), which are discussed in detail in section 3, they are qualitatively very similar, and the explanations for the behavior of ∥ S p 2 is essentially the same in both theories.
The effect of λ Re in the PT is qualitatively similar to that in the ZT; that is, it captures the classical effect of λ Re on the inertial subrange through its effect on S f 2 , τ TR and τ R .

Gustavsson et al theory
The Gustavsson et al theory [34] (hereafter GT) describes the relative velocities of inertial particles in isotropic turbulence, with particular attention given to the role played by the formation of caustics in the velocity distribution of the particles. In a similar manner to the ZT a transport equation for the PDF ( ) p t r w , , is considered, however an important difference is that in the GT, Δu is approximated as delta correlated in time. Consequently, the PDF equation they derive is a classical Fokker-Planck equation. From this they construct approximate analytical expressions for the moments of the PDF in the dissipation range of the form where the subscript N refers to the order of the moment, D 2 is the phase-space correlation dimension, and d is the physical space dimension (i.e. d = 3 for a 3D flow). The first term on the rhs of (16) describes the smooth part of the relative velocity moment and the second term describes the caustic contribution, i.e. that associated with the non-local particle dynamics. The theory is unclosed since f St ( ) The second-order structure function for the particles is found from (16)  is the configuration space correlation dimension, which then gives  It is interesting to note that their caustic contribution (second term on the rhs of (17)) is directly related to the spatial clustering via d St ( ) 2 , In the ZT equation for S p 2 (i.e. (5)) the non-local contribution (which gives rise to caustics) is also related to the spatial clustering through ρ. In general the GT predicts a bi-fractal behavior for the longitudinal relative velocities ∝ 2 and we note that a similar behavior has recently been derived for the behavior of the wall-normal particle velocities in wall-bounded turbulent flows [35].
For ) and finite relative velocities at r = 0. This is in contrast to the ZT.
Despite this disagreement, there are important similarities among the three theories. In particular, they all predict a scale transition from a smooth to a non-smooth region as r is decreased. Furthermore, in [10] it was shown using DNS data that the GT prediction that the structure function power law saturation exponent should be equal to − d d 2 for small r and large St is quite accurate for ≳ St 1. We shall return in section 4 to consider in more detail some of the scaling predictions of the GT and how they relate to the ZT.

Comparing theory with DNS
In this section, we compare the predictions for ∥ S r ( ) ). The ZT and PT solutions are generated by solving their equations given in [20,25] respectively. For details on the DNS see Ray and Collins [36].
The results are shown in figure 4. Before proceeding with the analysis, we note the striking qualitative similarity between these results for particle pairs in homogeneous, isotropic turbulence and those obtained for inertial particles in a turbulent boundary layer [28]. In particular, the near-wall accumulation and enhanced deposition velocity of particles in a boundary layer are somewhat analogous to the accumulation of pairs and their enhanced relative velocity in isotropic turbulence.
The results show that for = St 0.7, 1 the PT provides better agreement with the DNS than the ZT, both qualitatively and quantitatively, whereas for St = 3 the ZT gives the better predictions (in the dissipation range). For ⩽ St 1 the PT captures the significant departures from the smooth scaling in the dissipation regime better than the ZT. Indeed, in contrast to the DNS and PT, for St = 0.7 the ZT prediction for ⊥ S p 2 is very close to ⊥ S f 2 . From these comparisons with the DNS it is clear there is room for improvement in both theories. Below we discuss some causes of the errors. The PT consistently under-predicts S p 2 in the dissipation regime. Recall that the feature giving rise to > S S p f 2 2 is the non-local contribution to the particle velocity dynamics; pairs arriving at r may have relative velocities larger than the local fluid velocity difference since they are carrying the energy they received at larger separations in their history. The non-local contribution to S p 2 is therefore related to the backward-in-time separation of the pair. Recall from section 2.2 that although it was stressed in the PT that the backward-in-time separation of the pairs is an important feature to properly describe S p 2 , they modeled φ based on results for forward-in-time dispersion from Bec et al [9] (although they used the backward-in-time value for the Richardson constant in the RO regime, see section 2.2). Backward-in-time dispersion of inertial particle pairs is likely to have a faster rate than forward-in-time dispersion; this is certainly true for fluid particles [32,33]. Should this be the case, then the introduction of a backward-in-time model for φ would increase the PT predictions for S p 2 , potentially improving the agreement with the DNS. In addition, the PT does not incorporate the exponential separation behavior that occurs for fluid particle dispersion in the dissipation range, and probably also occurs for low St particles. Since the exponential separation is based on the fact that in the dissipation range Δ ∝ u r, and since this linear behavior is known to persist up to η ≈10 [23] then we would expect this to affect S p 2 up to a few Kolmogorov lengths, and its omission could explain some of the remaining discrepancy between the PT and DNS in the dissipation range (at least for lower Stokes numbers). We suspect that errors in the specification of φ are a significant source of error in the PT; this suggests a potential direction for improving the PT would be to introduce a backward-in-time model for φ into the PT.
Comparing the ZT predictions with DNS in figure 4 reveals that there are significant errors in the ZT predictions for St = 0.7 and 1.0, but reasonably good agreement for St = 3 (at least in the dissipation range; for η ≳ r 10 the ZT significantly overpredicts the DNS). There are several possible sources of error in the ZT discussed below.
The first potential source of error in the ZT is the QGA for S p 4 in the closure of S p 3 . Since the PDFs for inertial particle relative velocities in turbulence are strongly non-Guassian with heavy tails at the small scales (e.g. [27,36]), a QGA for S p 4 is likely to lead to errors in S p 2 in (4). We can use the DNS data to assess the implications of the QGA, under which and from (9) We will test (18) and (19) by using the DNS data for ∥ S p 2 to evaluate the right-hand sides, and then compare the predictions with the DNS data for ∥ S p 4 and ∥ S p 3 . The results are shown in figure 5. Figures 5(a), (c) and (e), which compare (18) with the DNS, show that (18) substantially under-predicts the DNS at the small scales, but becomes increasingly accurate with increasing r. This reflects the fact that ∥ p r w ( , ) is highly non-Gaussian at small scales [27,36], but tends toward a Gaussian distribution as → r L I . Figures 5(b), (d) and (f), which evaluate (19), show that (19) too has significant errors at small separations, for related reasons. Note that for St = 0.7 and 1, (19) under-predicts the DNS, whereas for St = 3, (19) over-predicts the DNS at the small scales. This is consistent with the behavior observed in figure 4, where the ZT under-predicts the DNS at the small scales for St = 0.7 and 1, but over-predicts it for St = 3.  (18)  N . An alternative to QGA would be to use the quadrature method of moments (QMOM, e.g. [37,38]) to close the moment equations in the ZT. This may preserve the correct scaling for ∥ S N p . There is work underway to explore this idea and, if successful, this will be the topic of a future paper. A second closely related error is the Gaussian assumption for Δu. For finite inertia particles, a memory of large values of Δu would generate large relative velocities as the particles approach each other. Therefore the neglect of the strongly non-Gaussian properties of Δu at the small scales could exacerbate the errors associated with the QGA.
In Part I we discussed the implications of an additional drift mechanism in the PDF equation-a drift mechanism that arises from the r dependence of the Δ ( ) t u r, statistics and their finite temporal correlation radius-that was incorrectly omitted in the ZT. In the moment equations, this drift arises explicitly in the RDF equation, but only implicitly in the equation for S p 2 through its effect upon the RDF. In Part I we showed using perturbation analysis that in the limits ≪ St 1 and η ≪ r , this drift is negligible, and using DNS data we inferred that it remains small even for ∼ St 1. In the inertial range this drift may not be zero, even in the limit ≪ St 1. However, we do not believe its omission is a source of significant errors in the inertial range shown in figure 4 because at these Stokes numbers, and at this modest Reynolds number, there is very little clustering occurring outside of the dissipation range.
The ZT (like the PT) approximates the time scales of the fluid velocity differences along inertial particle pair trajectories by those along fluid particle pair trajectories, another potential source of error. Modeling the effect of inertia on these time scales is non-trivial because of the way the particle pair is swept by the scales of the turbulence. At the large scales, the Eulerian time scale of Δu is larger than the Lagrangian time scale and so inertia would tend to increase the time scale of Δu experienced by the inertial particles. On the other hand, at the small scales, the Lagrangian time scales of Δu are associated with the Lagrangian time scales of the rate-ofstrain and rate-of-rotation tensors of the turbulence; strain and rotation are small scale quantities but the single-particle trajectories along which they are evaluated are strongly influenced by the larger scales of the flow. Therefore understanding the influence of St on the time scales of Δu in the dissipation range is non-trivial and will need to be investigated in future work.
A related error in the ZT results from the use of local closures for λ and μ even though they are inherently non-local due to the r dependence of the Δ t u r ( , ) statistics (see Part I for their unclosed definitions). Local closures for these dispersion tensors can be in considerable error [39]; for example, in Part I we introduced a non-local correction to the local closure for the diffusion coefficient λ in the RDF equation, and showed its importance for obtaining accurate predictions of the RDF. However in contrast to the RDF, the non-local corrections have a relatively weak effect on the ZT predictions for S p 2 . The reason is that the equation for S p 2 is influenced by μ much more than λ, and μ, unlike λ, is only weakly non-local [39].
It is clear from the above results that while the theories are able to describe the qualitative features of S p 2 , there is still some way to go before they will be capable of making quantitatively accurate predictions.
Finally  versus St 1 , as in the limit → r 0 this would yield a straight line according to (20); however, for finite η r the result is the curve shown in figure 6. In addition to the DNS, we show the predictions of the ZT, the PT and (21) (referred to in the plot as 'WM' for Wilkinson and Mehlig) with coefficients A = 1.2, = η B u 2.2 taken from [11].
Overall the PT provides the best prediction; it captures the qualitative trend of the DNS better than the ZT or (21), although as discussed earlier it tends to under-predict the DNS data. For the larger St numbers, one advantage the PT and ZT have over the WM is that the latter was derived in [40] under the assumption that Δ ∝ u r. While this smooth approximation is valid in the dissipation range, it is not valid in the inertial range, and particles with larger Stokes numbers carry some memory of their path histories from the non-smooth inertial range when their current separation is in the dissipation range. The PT and the ZT contain models for the inertial range and so, in principle, capture this effect, whereas the WM does not. We also note that (21) was derived making a white-in-time approximation for Δu p . In [41], the white-noise assumption is relaxed to take into account a small, but finite time correlation for Δu p . This relaxed assumption yielded The ZT does not accurately predict the DNS data in figure 6 over the entire Stokes number range. In particular, the caustic contribution to ⎡ ⎣ ⎤ ⎦ ∥ ( ) S r ln p 2 is not significant until ≳ St 1 and then it grows more rapidly than the DNS. The reason for this is not entirely clear. As discussed earlier, the non-local contribution in the ZT that gives rise to the caustic depends upon S p 3 and this is not correctly predicted because of the QGA. In [25], the authors note that built into the ZT is a mean square separation behavior that yields a Richardson constant g ≈ 3, approximately three times larger than the backward Richardson constant found in experiment [33] and used in the PT. Since only larger St particles at η ≲ r will be affected by their backward-in-time separation behavior at times in the past when they were in the inertial regime, this may in part explain why the ZT prediction for ⎡ ⎣ ⎤ ⎦ ∥ ( ) S r ln p 2 grows too fast with St for ≳ St 1. This probably also partially explains some of the overpredictions shown in figure 4 from the ZT for S p 2 in the inertial range.

Relationship between S p 2 and the RDF
We now wish to consider the impact the relative velocity structure functions have on the RDF, and in particular how caustics affect the RDF. The ZT gives a direct relation between S p 2 and the RDF through (4). When cast into spherical coordinates, and assuming isotropy, the solution to (4) can be written as

( )
It is worth highlighting that the scale break, which is a point of difference between the GT and the ZT, does not affect the scaling of density-weighted moments of the PDF. This is an important point since in some instances (for example the collision kernel), it is the density weighted moments of the PDF that are of interest. We define the density weighted moments of the longitudinal particle relative velocity PDF as where the r 2 factor comes from integration over the spherical angles (which is performed in the GT definition of the moments). Now for N = 2, m 2 scales like r 4  , m 2 scales like r 2 because the product ∥ gS p 2 is independent of r in this region. The interesting question is whether the scale break predicted by the ZT affects this behavior at small r. Figure 9 shows a plot of the ZT prediction for m 2 at = St 1.5, 2.5. As you can see, the m 2 predicted by the ZT scales like r 2 for η ≲ r for these values of St. Figure 8 shows that for St = 1.5 the scale break occurs at η ∼ − r 10 1 , however figure 9 (a) clearly shows that this has no effect upon m 2 , and the same is true of the St = 2.5 case. Thus, the scale break does not affect the behavior of m 2 , and consequently the GT and ZT predictions for m 2 are in qualitative agreement.

Conclusion
In this study, we compared three theoretical predictions for the relative velocity statistics of inertial particles embedded in an isotropic turbulent flow. The focus of the three theories is on the second-order structure function (both along the line of centers and perpendicular to the line of centers) as a function of the inertia of the particles as measured by the particle Stokes number based on the Kolmogorov time scale. For ≪ St 1, the theories predict smooth behavior for the structure functions; however, at larger Stokes numbers the theories display a non-smooth behavior in the dissipation range of the turbulence, which is related to the formation of caustics. The physical argument for the formation of these caustics is similar in the three theories; however, their quantitative predictions differ significantly. Note that the GT does not specify all of the coefficients of the theory, and so only qualitative information can be obtained.
Comparisons with DNS data show that while the ZT and PT are in qualitative agreement with each other and the DNS, the PT provides a more accurate prediction of the DNS. However, even for the PT the errors are still significant. One source of error in the PT is the use of forward-in-time (rather than backward-in-time) dispersion closures. It is known that backwardin-time dispersion is more rapid, and hence this may explain the under-predictions by the PT. Refining the approximations made in the PT is a challenge, however, we think that this could be done by incorporating a better understanding of the backward-in-time dispersion of inertial particle pairs.
Finally, we showed how the formation of caustics might impact the RDF by considering the general expression for the RDF from the ZT. Note that this expression was shown to be accurate in Part I of this paper. We find that the non-local contribution to the particle velocity dynamics can induce a scale break in the RDF, causing the power-law behavior to transition to a constant plateau at very small separations. For sufficiently small Stokes numbers, this transition occurs at separations below the particle size, making it irrelevant from a practical sense, as other neglected physics (e.g. hydrodynamic interactions or particle collisions) would likely become important before the plateau region is reached. However, for ∼ St 1, the transition could be important and would tend to reduce the RDF relative to the pure power-law form. Unfortunately, it is difficult to test this prediction with DNS, as it requires accurate statistics at separations well below the Kolmogorov length scale, and this translates to the requirement of an extremely large population of particles, which is challenging computationally. However, at still larger Stokes numbers this transition occurs in the inertial range, well within the reach of DNS. Results in the literature for such particles confirm qualitatively the ZT prediction of how caustics modify the scaling of the RDF.
This predicted scale break in the RDF coincides with the second-order structure function transitioning from a power law to a constant at small r, something predicted by the ZT and PT but not by the GT. We have argued that the reason for this discrepancy is that the GT makes the a priori assumption that the RDF is a power law in the limit → r 0, ruling out the possibility of a scale break in the particle velocity structure functions. This point of difference between the theories is something that needs to be resolved in future work.