Relative velocities in bidisperse turbulent suspensions

TL;DR

This paper studies how the relative velocities of particles of different sizes in turbulence are distributed, revealing power-law tails and their dependence on particle inertia, with implications for understanding particle interactions.

Contribution

It introduces a statistical model for bidisperse turbulent suspensions and analyzes the distribution of relative velocities, highlighting non-analytic behavior of power-law exponents.

Findings

Power-law tails occur when Stokes numbers are similar.

Exponent depends non-analytically on mean Stokes number.

Tail dominance persists even with larger Stokes-number differences.

Abstract

We investigate the distribution of relative velocities between small heavy particles of different sizes in turbulence by analysing a statistical model for bidisperse turbulent suspensions, containing particles with two different Stokes numbers. This number, ${\rm St}$, is a measure of particle inertia which in turn depends on particle size. When the Stokes numbers are similar, the distribution exhibits power-law tails, just as in the case of equal ${\rm St}$. The power-law exponent is a non-analytic function of the mean Stokes number $\overline{\rm St}$, so that the exponent cannot be calculated in perturbation theory around the advective limit. When the Stokes-number difference is larger, the power law disappears, but the tails of the distribution still dominate the relative-velocity moments, if $\overline{\rm St}$ is large enough.

Equations

$$\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{x}$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\Delta x=\Delta v\text{mbox}and\frac{\text{d}}{\text{d}t}\left[\begin{array}{c}\Delta v\\ 2\overline{v}\end{array}\right]=\left[\begin{array}{cc}1& \theta \\ \theta & 1\end{array}\right]\left[\begin{array}{c}A(t)\Delta x-\Delta v\\ B(t)-2\overline{v}\end{array}\right],$$

$$g_{\mu }(\Delta x)=|\Delta x{|}^{\mu -1}\text{mbox}and\frac{\mathrm{d}}{\mathrm{d}z}(z+z^2+{\overline{\epsilon }}^2\frac{\mathrm{d}}{\mathrm{d}z})Z_{\mu }=\mu z{Z}_{\mu }$$

$$Z_{{\mu }_{\mathrm{c}}}(z)=Z^{(0)}(z)+\delta \mu Z^{(1)}(z)+\delta {\mu }^2{Z}^{(2)}(z)+\dots .$$

$$Z_{-1}=C_1(z+1){\mathrm{e}}^{-U(z)}+C_2(z+1){\int }_{-\infty }^z\text{d}t\frac{{\mathrm{e}}^{U(t)-U(z)}}{(t+1{)}^2}.$$

$$Z^{(1)}=\frac{z+1}{{\overline{\epsilon }}^2}{\int }_{-\infty }^z\text{d}t\frac{{\mathrm{e}}^{U(t)-U(z)}}{(t+1{)}^2}{\int }_{-1}^t\text{d}{t}^{\prime }{t}^{\prime }(t^{\prime }+1)Z_{-1}(t^{\prime }).$$

$${\mu }_{\mathrm{c}}=-1+\delta \mu \text{mbox}with\delta \mu ={\pi }^{-1}\exp [-1/(6{\overline{\epsilon }}^2].$$

$$\frac{\mathrm{d}}{\mathrm{d}t}z=-\gamma z-z^2+\gamma A(x,t).$$

Full text

The dynamics of small heavy particles in turbulence plays a crucial role in many scientific problems and technological applications. Any model of the particle dynamics must refer to the turbulence the particles experience. This is a challenge for realistic modelling of such systems, for their direct numerical simulation (DNS), and for experiments. Novel particle-tracking techniques and improved DNS algorithms have made it possible to uncover striking phenomena in turbulent aerosols. For example, heavy particles tend to avoid the vortices of the turbulent fluid, and they form small-scale fractal patterns. Nearby particles can have very high relative velocities, an effect caused by ‘caustic’ singularities in the particle dynamics. The analysis of statistical models has led to substantial progress in explaining these phenomena, reviewed in Ref. Gustavsson and Mehlig (2016). This analysis has offered fundamental insights about how caustics shape the distribution of relative velocities of nearby particles Sundaram and Collins (1997); Falkovich et al. (2002); Wilkinson and Mehlig (2005); Wilkinson et al. (2006); Bec et al. (2010); Gustavsson and Mehlig (2011); Salazar and Collins (2012); Bewley et al. (2013); Perrin and Jonker (2015).

These results apply only to ‘monodisperse’ suspensions of identical particles. The question is therefore how particles of different sizes cluster and move relative to each other. Spatial clustering of such ‘bidisperse’ suspensions was analysed in Refs. Chun et al. (2005); Bec et al. (2005). It was found that the particles cluster onto two distinct attractors, and that the fractal distribution of separations between differently-sized particles is cut off at a small spatial scale, $r_{\rm c}$, that can be much larger than the particle size.

In this Letter we analyse the distribution of relative velocities of particles with different sizes. Our results are important for the physics of turbulent aerosols, because the distribution of relative velocities determines the speeds of colliding particles, their collision rate, and collision outcomes Wilkinson et al. (2008); Windmark et al. (2012); Gustavsson and Mehlig (2014).

In a dilute suspension of small, heavy, spherical particles the dynamics of a single particle is approximately given by Stokes law with constant $\gamma\equiv 9\rho_{f}\nu/(2a^{2}\rho_{p})$:

[TABLE]

Here $\boldsymbol{x}$ and $\boldsymbol{v}$ are particle position and velocity, $a$ is the particle size, $\nu$ is the fluid viscosity, and $\rho_{f}$ and $\rho_{p}$ are fluid and particle densities. We model the turbulent fluid velocities by a random Gaussian velocity field $\boldsymbol{u}(\boldsymbol{x},t)$ with zero mean, correlation time $\tau$, correlation length $\eta$, and typical speed $u_{0}$, representing the universal small spatial scales of turbulence Schumacher et al. (2014), neglecting intermittent, non-Gaussian features Gustavsson and Mehlig (2016). Eq. (1) assumes that the particle and shear Reynolds numbers are small. Gravitational settling Gustavsson et al. (2014); Bec et al. (2014); Ireland et al. (2016); Mathai et al. (2016); Parishani et al. (2015) is disregarded, this is valid when the turbulence is intense enough. The dimensionless parameters of the model are the Stokes number $\text{St}\equiv 1/(\gamma\tau)$, a measure of particle inertia, and the Kubo number, $\text{Ku}\equiv u_{0}\tau/\eta$, measuring the persistence of the flow.

Statistics of relative velocities and separations. Fig. 1a shows results of statistical-model simulations for $\text{Ku}=1$ in two spatial dimensions. Panel a shows the distribution $\varrho(v_{r},r)$ of relative radial velocity $v_{r}$ and separation $r$ between particles with slightly different Stokes numbers. For fixed small $v_{r}$, the distribution exhibits a power law as a function of $r$ for large values of $r$, and becomes uniform for small values of $r$. This cross over defines the cut-off scale $r_{\rm c}$, observed and discussed in Refs. Chun et al. (2005); Bec et al. (2005).

Fig. 1a shows that there is a second cut-off scale, $v_{\rm c}$, that distinguishes between a power law as a function of $v_{r}$ for large $v_{r}$, and a plateau for small $v_{r}$. Our numerical results demonstrate that this new scale $v_{\rm c}$ does not depend on $r$ (Fig. 1a). This is in marked difference to the monodisperse case, where the distribution has no plateau in the limit of $r\to 0$. Panel b shows that $v_{\rm c}$ depends linearly on the parameter $\theta\equiv|\text{St}_{1}-\text{St}_{2}|/(\text{St}_{1}+\text{St}_{2})$ measuring the difference between the Stokes numbers. Panel c shows distributions for $r\ll r_{\rm c}$ as functions of $v_{r}$, for different values of $\theta$. The dashed line is the power law ${|v_{r}|}^{\mu_{\rm c}-d-1}$, where $d=2$ is the spatial dimension. Panel d demonstrates that the exponent $\mu_{\rm c}$ is approximately equal to ${\rm min}(D_{2},d+1)$ where $D_{2}(\overline{\text{St}})$ is the phase-space correlation dimension of a suspension of identical particles with Stokes number $\overline{\text{St}}\equiv(\overline{\gamma}\tau)^{-1}$ and $\overline{\gamma}=\tfrac{1}{2}(\gamma_{1}+\gamma_{2})$.

The shape of the distribution shown in Fig. 1c has a strong effect on the moments of $v_{r}$. For bidisperse suspensions, the uniform part of $\varrho(v_{r},r)$ ($r<r_{\rm c}$ and $|v_{r}|<v_{\rm c}$) yields an $r^{d-1}$-contribution to the moments $\langle|v_{r}|^{p}\rangle$, just as the caustic-contribution for identical particles Gustavsson et al. (2012), caused by uniformly distributed, uncorrelated particles with large relative velocities originating from caustic singularities Falkovich et al. (2002); Wilkinson and Mehlig (2005). For small mean Stokes number, $\overline{\text{St}}\ll 1$, when caustics are rare, the dominant contribution to the moments of $v_{r}$ at small $r$ comes from the uniform part of the distribution. But for $\overline{\text{St}}\sim 1$, caustics are abundant and for $p\geq 1$ the moments $\langle|v_{r}|^{p}\rangle$ are dominated by the tails of the distribution.

The uniform part of the bidisperse distribution below $r_{\rm c}$ and $v_{\rm c}$ affects the tails indirectly. Comparing the normalised mono- and bidisperse distributions we see that the tails for $\theta\neq 0$ must lie above those for $\theta=0$. This means that the moments $\langle|v_{r}|^{p}\rangle$ must have a minimum for $\theta=0$ (at fixed $\overline{\text{St}}$ and fixed small separation). This is consistent with the findings of Refs. Völk et al. (1980); Mizuno et al. (1988); Markiewicz et al. (1991); Pan et al. (2014a); James and Ray (2016).

Analysis of the white-noise limit. We now show that all of the above observations can be explained qualitatively by analysing a one-dimensional white-noise model. We linearise Eq. (1), $\Delta u(x,t)\sim A(t)\Delta x$, to obtain equations for small separations $\Delta x\equiv x_{2}-x_{1}$. The spatial dependence in $A(x,t)\equiv$ $\partial_{x}u(x,t)$ is neglected here, this disregards preferential sampling which is not important in the white-noise limit Gustavsson and Mehlig (2016). The gradient $A(t)$ has zero mean and correlation function $\langle A(t_{1})A(t_{2})\rangle=3(u_{0}/\eta)^{2}\exp(-|t_{2}-t_{1}|/\tau)$. For $\theta\neq 0$ there is an additional stochastic driving, $u(x_{1},t)+u(x_{2},t)\approx 2u(\overline{x},t)$ where $\overline{x}\equiv\tfrac{1}{2}(x_{1}+x_{2})$. Neglecting the spatial dependence we write $B(t)\equiv 2u(\overline{x},t)$. The two noise terms $A(t)$ and $B(t)$ are uncorrelated, and $\langle B(t_{1})B(t_{2})\rangle=4u_{0}^{2}\exp(-{|t_{2}-t_{1}|}/{\tau})$. We use the de-dimensionalisation $\tilde{t}\equiv t\overline{\gamma}\,,\tilde{x}\equiv x/\eta\,,\tilde{v}\equiv v/(\eta\overline{\gamma})\,,\tilde{u}\equiv u/(\eta\overline{\gamma})$. Dropping the tildes we find:

[TABLE]

where $\Delta v\equiv v_{2}-v_{1}$ is the relative velocity. The white-noise limit is taken by letting $\overline{\text{St}}\to\infty$ and $\text{Ku}\to 0$ such that $\overline{\varepsilon}^{2}\equiv 3\text{Ku}^{2}\,\overline{\text{St}}$ stays finite. In this limit we obtain $\langle A(t_{1})A(t_{2})\rangle=\tfrac{3}{4}\langle B(t_{1})B(t_{2})\rangle=2{\overline{\varepsilon}}^{2}\delta(t_{2}-t_{1})$. The white-noise problem has three coupled dynamical variables, $\Delta x$ and $\Delta v$ as in the monodisperse case, and the mean velocity $\overline{v}\equiv{\tfrac{1}{2}}(v_{1}+v_{2})$. We write $\varrho(\Delta v,\Delta x)\equiv\int{\rm d}\overline{v}\,P(\Delta v,\Delta x,\bar{v})$, where $P(\Delta v,\Delta x,\bar{v})$ is the steady-state distribution of $\Delta v,\Delta x$ and $\bar{v}$.

Cross-over scale $v_{\rm c}$. For small values of $\theta$, the mean velocity $\overline{v}$ is approximately governed by an Ornstein-Uhlenbeck process driven by $B$. In the equation for $\Delta v$ we identify three competing terms: the two noise terms $V_{1}\equiv A\Delta x$ and $V_{2}\equiv\theta(B-2\overline{v})$, and the damping $-\Delta v$.

For small relative velocities ($|\Delta v|\to 0$), the damping term is negligible. In this limit, the nature of the dynamics depends on whether $|\Delta x|\ll\theta$ or $|\Delta x|\gg\theta$. In the first case $V_{2}$ dominates over $V_{1}$, in the second case $V_{1}$ dominates. We find the spatial cut off $r_{\rm c}$ by estimating the value of $\Delta x$ where $V_{1}$ and $V_{2}$ are of the same order. This yields $r_{\rm c}\propto\sqrt{\Phi_{V_{2}}/\Phi_{A}}\propto\theta$, where $\Phi_{X}$ is the integrated correlation function of $X$. Now consider small separations ($|\Delta x|\ll r_{\rm c}$) and finite $\Delta v$. In this case $V_{1}$ is negligible. Now we must consider whether $|\Delta v|\ll\theta$ or $|\Delta v|\gg\theta$. In the first case the noise term $V_{2}$ dominates, so that the dynamics is diffusive. When $|\Delta v|\gg\theta$, on the other hand, the dynamics is deterministic, and the noise plays no role. Between these two regimes we find the cut-off scale $v_{c}\propto\sqrt{\Phi_{V_{2}}}\propto\theta$. These considerations explain the linear $\theta$-dependencies of the cutoff-scales $r_{\rm c}$ and $v_{\rm c}$, Figs. 1c and Fig. 2a. Finally, consider the case where $V_{2}$ is negligible ($|\Delta x|\gg r_{\rm c}$ or $|\Delta v|\gg v_{\rm c}$). In this limit the dynamics is essentially that of monodisperse suspensions, but with mean Stokes number $\overline{\text{St}}$. For $|\Delta x|>r_{\rm c}$ or $|\Delta v|>v_{\rm c}$ the distribution is expected to show the same power-law tails as the monodisperse case (Fig. 1c).

Power-law tails. We now discuss how the power-law tails observed in bidisperse suspensions for small $\theta$ are related to the power laws observed in the monodisperse case. For $\theta\ll 1$ we decompose the joint distribution $P(\Delta v,\Delta x,\bar{v})=\varrho_{0}(\Delta v,\Delta x)p(\overline{v})+\delta P(\Delta v,\Delta x,\overline{v})$, where $\varrho_{0}(\Delta v,\Delta x)$ equals the distribution for monodisperse particles ($\theta=0)$. It follows from the analysis above that the relative dynamics is diffusive for $|\Delta x|\ll r_{\rm c}$ and $|\Delta v|\ll v_{\rm c}$. In this case $\delta P(\Delta v,\Delta x,\bar{v})$ cannot be neglected. For $|\Delta x|\gg r_{\rm c}$ or $|\Delta v|\gg v_{\rm c}$, by contrast, the dynamics is that of the monodisperse suspension and $\varrho_{0}(\Delta v,\Delta x)$ dominates over $\delta P(\Delta v,\Delta x,\bar{v})$. To study the tails it thus suffices to consider $\varrho_{0}(\Delta v,\Delta x)$. We follow Ref. Gustavsson and Mehlig (2011) and make the separation ansatz: $\varrho_{0}(\Delta x,z)=\sum_{\mu}A_{\mu}g_{\mu}(\Delta x)Z_{\mu}(z)$ with $z\equiv\Delta v/\Delta x$ . Inserting this ansatz into the Fokker-Planck equation corresponding to Eq. (2) we obtain

[TABLE]

to lowest order in $\theta$. This is the equation for a suspension of identical particles Gustavsson and Mehlig (2011) with Stokes number $\overline{\text{St}}$. In a slightly different form, Eq. (3) was used to model the effect of wall collisions in turbulent suspensions Belan et al. (2014).

Particle-exchange symmetry Gustavsson and Mehlig (2011) requires that $Z_{\mu}(z)$ is symmetric for large $|z|$, $\lim_{z\to\infty}Z_{\mu}(z)/Z_{\mu}(-z)=1$. Numerical analysis shows that Eq. (3) exhibits a discrete set of such solutions Gustavsson and Mehlig (2014). For small $\overline{\varepsilon}$ only two values are allowed: $\mu=0$ and $\mu=\mu_{\rm c}(\overline{\varepsilon})$. At large values of $|z|$ one finds that $Z_{\mu_{\rm c}}(z)\sim|z|^{\mu_{\rm c}-2}$. In the monodisperse case, $\mu_{\rm c}({\overline{\varepsilon}})$ equals $D_{2}$ for $\overline{\varepsilon}>\overline{\varepsilon}_{\rm c}$. For $\overline{\varepsilon}<\overline{\varepsilon}_{\rm c}$, by contrast, particle paths coalesce exponentially Wilkinson and Mehlig (2003), there is no power law steady-state, and $\mu_{\rm c}({\overline{\varepsilon}})$ becomes negative and loses its meaning as correlation dimension. The distribution $g_{\mu_{\rm c}}(\Delta x)=|\Delta x|^{\mu_{\rm c}-1}$ is not normalisable for $\mu_{\rm c}\leq 0$.

In the bidisperse case, this divergence is regularised since diffusion dominates for $|\Delta x|\ll r_{\rm c}$ and $|\Delta v|\ll v_{\rm c}$. More precisely, $\delta P(\Delta v,\Delta x,\bar{v})$ regularises $P(\Delta v,\Delta x,\overline{v})$. Thus a power-law steady-state distribution is obtained, even for $\overline{\varepsilon}<\overline{\varepsilon}_{\rm c}$. This is analogous to the effect of small-scale diffusion that regularises a power-law steady state distribution of separations Wilkinson et al. (2015).

The full solution of (3) consists of a sum of two terms, $A_{0}\Delta x^{-1}Z_{0}(z)+A_{\mu_{\rm c}}|\Delta x|^{\mu_{\rm c}-1}Z_{\mu_{\rm c}}(z)$. As was discussed in Gustavsson and Mehlig (2011), the separation ansatz for $\varrho(\Delta x,z)$ is strictly valid only for $|\Delta x|\ll 1$. The $A_{0}$-term is therefore subdominant for $\mu_{\rm c}<0$. For $\mu_{\rm c}>0$ we expect $A_{0}\sim\theta\ll A_{\mu_{\rm c}}$ since $A_{0}$ must vanish in order for the distribution to be normalisable over $\Delta x=0$ as $\theta\to 0$, consistent with the numerical small-$\theta$ results in Fig. 2a. Therefore we disregard the $\mu=0$-term for all values of $\overline{\varepsilon}$.

Fig. 2b shows that for $\theta\ll 1$ the exponents of both the $\Delta x$- and $\Delta v$- distributions (symbols) follow the numerically calculated $\mu_{c}(\bar{\varepsilon})$ very precisely (solid line). The asymptote $Z_{\mu_{\rm c}}(z)\sim|z|^{\mu_{\rm c}-2}$ for large $|z|$ gives rise to the power-law tails in the relative-velocity distribution of the form $|\Delta v|^{\mu_{\rm c}-2}$. In $d$ spatial dimensions a similar argument yields tails on the form $v_{r}^{\mu_{\rm c}-d-1}$. This explains the power laws observed in Figs. 1c and 2a.

In summary the white-noise analysis qualitatively explains not only the existence of the plateau in the relative-velocity distribution and the corresponding cutoff scales $r_{c}$ and $v_{c}$. It also yields the linear dependence of these scales upon $\theta$, and explains the power-law tails. Thus, the one-dimensional white-noise analysis qualitatively explains all the observed features in Fig. 1.

Failure of perturbation theory in $\overline{\varepsilon}$. Fig. 2b demonstrates that the exponent $\mu_{\rm c}$ depends very sensitively on ${\overline{\varepsilon}}$, for small values of ${\overline{\varepsilon}}$. Here we show that this dependence is non-analytic. This is important, because it means that perturbation theory in ${\overline{\varepsilon}}$ fails. This is, in turn, a likely reason Gustavsson et al. (2015) why Borel resummation of perturbation theory does not yield accurate results for the correlation dimension $D_{2}$ (Fig. 1 in Ref. Wilkinson et al. (2010)). There is to date no theory explaining these observations.

The exponent $\mu_{\rm c}$ appears as a generalised eigenvalue in Eq. (3), but exact solutions to Eq. (3) are known only for $\mu=0$ Wilkinson and Mehlig (2003), $\mu=-1$ Derevyanko et al. (2007), and in the large $\overline{\varepsilon}$-limit Schomerus and Titov (2002). The physical solution at small values of ${\overline{\varepsilon}}$ corresponds to the eigenvalue $\mu_{\rm c}({\overline{\varepsilon}})$. The solution for $\mu=-1$ is unphysical since it does not obey particle-exchange symmetry.

Numerical analysis reveals that $\mu_{\rm c}({\overline{\varepsilon}})\to-1$ as ${\overline{\varepsilon}}\to 0$. We write $\mu_{\rm c}\!=\!-1\!+\!\delta\mu$, where $\delta\mu\geq 0$ is an ‘eigenvalue splitting’ that vanishes as $\overline{\varepsilon}\to 0$, and attempt to find the physical solution by perturbation theory in $\delta\mu$:

[TABLE]

Substituting (4) into (3) we obtain a hierarchy of differential equations for $Z^{(n)}$ (Supplemental Material Supplemental material to this article (2016)). The lowest-order solution is the unphysical one, $Z^{(0)}=Z_{-1}$. Its general form is known Derevyanko et al. (2007):

[TABLE]

Here $U(z)\equiv(\tfrac{1}{3}z^{3}+\tfrac{1}{2}z^{2})/{\overline{\varepsilon}}^{2}$, and $C_{1}$ and $C_{2}$ are constants. The first-order solution reads:

[TABLE]

The integrand involves $Z_{-1}$, and we must take $C_{1}=0$ in Eq. (5) for $Z_{-1}$ since the first term in (5) is not normalisable. We further choose $C_{2}=-{\overline{\varepsilon}}^{2}$ because it normalises the large-$z$ tails of $Z_{-1}$ to $Z_{-1}\sim-z^{-3}$. To evaluate (6) we use the fact that $\delta\mu\to 0$ as $\overline{\varepsilon}\to 0$ and apply a WKB approximation Bender and Orszag (1978) of $Z_{-1}$ that becomes exact in the limit $\overline{\varepsilon}\to 0$. The equation for $Z_{-1}$ has two ‘turning points’ at $z=-1,0$ where the WKB approximation breaks down. We find locally exact solutions and use these to match the WKB solutions across the turning points. Details of this asymptotic-matching method are given in the Supplemental Material Supplemental material to this article (2016). In this way we find $Z^{(1)}\sim|z|^{-3}\log|z|$ for large negative values of $z$, and $Z^{(1)}\sim|z|^{-3}[2\pi{\rm e}^{1/(6{\overline{\varepsilon}}^{2})}-\log|z|]$ for large positive $z$. Imposing particle-exchange symmetry yields

[TABLE]

This result is shown as a dashed line in Fig. 2b. It is in good agreement with the numerical data for small ${\overline{\varepsilon}}$.

We emphasise that the non-analytic dependence (7) cannot be obtained by the ${\overline{\varepsilon}}$-perturbation theory described in Refs. Mehlig and Wilkinson (2004); Duncan et al. (2005); Gustavsson and Mehlig (2016). That approach starts by rescaling $z\to\overline{\varepsilon}z$, and solves the resulting equation $\tfrac{\rm d}{{\rm d}z}(z+\overline{\varepsilon}z^{2}+\tfrac{\rm d}{{\rm d}z})Z_{\mu}\!=\overline{\varepsilon}\mu zZ_{\mu}$ perturbatively in ${\overline{\varepsilon}}$. The solution remains essentially Gaussian. Even for very small ${\overline{\varepsilon}}$, however, the $z$-dynamics can escape from $z\approx 0$ to $-\infty$ by forming caustics Gustavsson and Mehlig (2016). This gives rise to power-law tails that are not described by the local theory. Higher orders in ${\overline{\varepsilon}}$ improve the accuracy only locally, inside an ${\overline{\varepsilon}}$-sized ‘boundary layer’ Hinch (1991) around $z=0$, but not in the tails. Consequently, perturbation theory fails to give the non-analytical dependence (7).

Our theory, by contrast, provides a uniform approximation valid for small $\overline{\varepsilon}$ and arbitrary $z$. It allows to compute the exponentially small eigenvalue splitting $\delta\mu$ because it captures the formation of caustics. For small ${\overline{\varepsilon}}$, the rate $J$ of caustic formation is $J\sim(2\pi)^{-1}{\rm e}^{-1/6{\overline{\varepsilon}}^{2}}$ Gustavsson and Mehlig (2016). It follows that $\mu_{\rm c}\sim 2J-1$ at small values of ${\overline{\varepsilon}}$.

We now show that caustics cause similar problems in small-St expansions Maxey (1987); Elperin et al. (1996); Balkovsky et al. (2001); Chun et al. (2005); Wilkinson et al. (2007); Gustavsson and Mehlig (2016). In one dimension, the equation for $z$ reads in dimensional units

[TABLE]

The Stokes number used in DNS of heavy particles in turbulence Bec et al. (2006, 2007) is $\kappa\equiv u_{0}/(\eta\gamma)$ Gustavsson and Mehlig (2016). Therefore we consider the small-$\kappa$ limit of (8), taking $A$ to be constant (in this limit the local fluid-velocity gradients are roughly constant during the time $\eta/u_{0}$ Gustavsson and Mehlig (2016)). As $\kappa\to 0$ we find that $z\sim A$, with distribution $P(A)=(6\pi u_{0}^{2}/\eta^{2})^{-1/2}\exp[-A^{2}/(6u_{0}^{2}/\eta^{2})]$, localised around zero. As in the white-noise case, this local approximation fails when caustics allow $z$ to escape to $-\infty$. Kramers-escape theory for coloured noise Tsironis and Grigolini (1988); Bray and McKane (1989) reveals that caustics occur in the persistent limit when $A<-\gamma/4$, with probability $p\!=\!\int_{-\infty}^{-\gamma/4}\!\!\text{d}A\,P(A)$. Evaluating this integral for $\kappa\ll 1$, one obtains a non-analytic dependence, $p\propto\exp[-1/(96\kappa^{2})]$, consistent with Ref. Gustavsson and Mehlig (2013). In DNS a similar non-analytic dependence is found Falkovich and Pumir (2007), albeit of a slightly different form because the turbulent velocity-gradients are non-Gaussian. We conclude that the perturbation theories in the white-noise limit and in the persistent limit fail for the same reason. Both expansions are valid only near $z\!=\!0$ and fail in the presence of caustics.

Conclusions. We analysed the distribution of relative velocities in turbulence between small, heavy particles with different Stokes numbers St. We demonstrated that the difference in St causes diffusive relative motion at small separations, giving rise to a plateau in the distribution for relative velocities smaller than a cut off $v_{\rm c}$. This is in qualitative agreement with DNS Pan et al. (2014b). Fig. 2 in Ref. Pan et al. (2014b) shows the $v_{r}$-distribution for a bidisperse suspension. The cut off $v_{\rm c}$ depends linearly on $\theta$, the parameter characterising the difference in Stokes numbers.

At small values of $\theta$ the distribution exhibits algebraic tails $|v_{r}|^{\mu_{\rm c}-d-1}$ as in the monodisperse case. The exponent $\mu_{\rm c}$ is determined by the phase-space correlation dimension of a monodisperse system with mean Stokes number $\overline{\text{St}}$. When $\overline{\text{St}}$ is $\mathcal{O}(1)$ or larger, the moments of relative velocities at small separations are dominated by the tails of the distribution. At small $\theta$ the tails are power laws, but as $\theta$ becomes larger, our statistical-model simulations show that the tails gradually disappear, consistent with $v_{\rm c}\propto\theta$. The DNS of Ref. Pan et al. (2014b) are for quite large $\theta$. It would be of interest to perform DNS for small values of $\theta$ to find the power-law tails we predict here.

Our analysis shows how sensitive the distribution is to polydispersity. This is important for experiments tracking the dynamics of $\mu$m-size particles in turbulence Bewley et al. (2013), where strict monodispersity is difficult to achieve.

We explained these observations using a one-dimensional statistical model. The difference in Stokes numbers regularises the distribution, so that the bidisperse model has a power-law steady-state at small $\overline{\varepsilon}$. This makes it possible to compute how the power-law exponent $\mu_{\rm c}$ depends on the inertia parameter $\overline{\varepsilon}$. We demonstrated that the dependence of $\mu_{\rm c}$ upon $\overline{\varepsilon}$ is not analytic.

Our theory explains why small-${\overline{\varepsilon}}$ expansions fail to give non-analytic contributions to physical quantitites in the white-noise limit. It is likely that non-analytic terms in white-noise expansions for the correlation dimension Gustavsson et al. (2015) and Lyapunov exponents Mehlig and Wilkinson (2004); Duncan et al. (2005) in two and three dimensions have similar origins, and we speculate that Eq. (7) contains only the first two terms in an infinite series of the form Dorigoni (2014) $\sum_{k=0}^{\infty}{\rm e}^{-k/(6{\overline{\varepsilon}}^{2})}\sum_{m=0}^{\infty}b_{m}^{(k)}{\overline{\varepsilon}}^{2m}$. There could also be logarithmic terms, $\log^{n}{\overline{\varepsilon}}^{2}$.

More generally we explained that small-St expansions for heavy particles in turbulence Maxey (1987); Elperin et al. (1996); Balkovsky et al. (2001); Chun et al. (2005); Wilkinson et al. (2007); Gustavsson and Mehlig (2016) suffer from similar problems when caustics occur. This indicates that matched asymptotic expansions (or similar methods) are required to explain the characteristic minimum Bec et al. (2007) of the correlation dimension as a function of St.

Our predictions can be directly tested by experiments or by DNS of turbulent bidisperse turbulent suspensions. We note, however, that our analysis pertains to the dynamics in the dissipative range. At higher Stokes numbers, when separations between particle pairs explore the inertial range Mehlig et al. (2007); Gustavsson et al. (2008), we expect corrections to the power-law exponents derived here.