1. Introduction
Collisions between bubbles and particles in a turbulent flow are of significant technological relevance. In particular, such collisions are essential to the flotation process, which is a widely used separation technique especially in the mining industry (Nguyen & Schulze Reference Nguyen and Schulze2004). In this process, after grinding, small ore fragments are fed into a big water-filled cell that is agitated via a rotor and into which air bubbles are injected. Collisions between the ore particles and the bubbles then form the base for the decisive test: valuable mineral particles attach to the bubbles due to their hydrophobic surface and consequently rise to the top where they can be skimmed off as a froth, whereas the hydrophilic waste rock particles remain in suspension and are eventually discharged as tailings. This technology is already applied at staggering scales (Nguyen & Schulze (Reference Nguyen and Schulze2004) estimated that a total of 2 billion tons of ore are treated annually). Given especially the relevance in the mining of copper, and in view of the strong push for electrification in response the the climate crisis, these numbers are likely to continue to rise in the future (Rogich & Matos Reference Rogich and Matos2008; World Bank Group 2017). The interest to understand the collision process better is driven by the demand for more reliable process modelling (Kostoglou, Karapantsios & Evgenidis Reference Kostoglou, Karapantsios and Evgenidis2020a) but also by the need for performance improvement. The latter is especially a concern for small particles with diameters smaller than 
$20\,\mathrm {\mu }{\rm m}$, where recovery is poor owing to their low collision rates (Nguyen, George & Jameson Reference Nguyen, George and Jameson2006; Miettinen, Ralston & Fornasiero Reference Miettinen, Ralston and Fornasiero2010).
For modelling purposes, the collision process is generally separated into two components (Pumir & Wilkinson Reference Pumir and Wilkinson2016): the ‘geometric collision rate’, which considers the collisions neglecting any interaction between the collision partners, and the ‘collision efficiency’, which quantifies how many of these collisions actually happen when taking the local modification of the flow field into account. The focus here is on the (ensemble-averaged) geometric collision rate between two species, ‘1’ and ‘2’, which when expressed per unit volume can be written as
\begin{equation} Z_{12} = \varGamma_{12} n_1 n_2, \end{equation}
where 
$\varGamma _{12}$ is the collision kernel, and 
$n_1$ and 
$n_2$ denote the respective number densities of the two species. The collision kernel measures the rate at which the separation vector between particle centres crosses the collision distance. For spherical particles with a collision distance 
$r_c = r_1+ r_2$ (where 
$r_1,r_2$ denote the particle radii), 
$\varGamma _{12}$ can be expressed as (Sundaram & Collins Reference Sundaram and Collins1997)
\begin{equation} \varGamma_{12} = 4{\rm \pi} r_c^2\,g(r_c)\,S_-(r_c), \end{equation}
where besides the surface area 
$4{\rm \pi} r_c^2$ of the collision sphere, the other factors are the radial distribution function (RDF) at collision distance 
$g(r_c)$, which describes variations in the local particle concentration, and the effective radial approach velocity at contact,
\begin{equation} S_-(r_c) ={-}\int_{-\infty}^{0} \Delta v_r\mathrm{p.d.f.}(\Delta v_r|r_c)\,\mathrm{d}(\Delta v_r). \end{equation}
Here, 
$\Delta v_r$ is the radial component of the relative velocity, which is positive when the pair separates, and 
$\mathrm {p.d.f.}(\Delta v_r|r_c)$ is the probability density function of 
$\Delta v_r$ conditioned on a pair with separation 
$r_c$.
The bidisperse collision kernel 
$\varGamma _{12}$ depends on a multitude of parameters characterising properties of the suspended particles and of the carrier flow. The discussion here is restricted to homogeneous isotropic turbulence for which the most relevant dependencies may be summarised in non-dimensional form as
\begin{equation} \frac{\varGamma_{12}}{r_c^3/\tau_\eta} = f\left(St_1,St_2,\frac{\rho_1}{\rho_f},\frac{\rho_2}{\rho_f},Re_\lambda, \frac{r_1}{\eta},\frac{r_2}{\eta},Fr,\ldots \right). \end{equation}Here, the Stokes number
\begin{equation} St_i = \frac{\tau_i}{\tau_\eta} = \frac{r_i^2(2\rho_i/\rho_f + 1)}{9\nu\tau_\eta}, \end{equation}
(
$i = 1,2$) characterises how well particles follow the flow by relating the particle response time 
$\tau _i = {r_i^2(2\rho _i/\rho _f + 1)}{/(9\nu )}$ to the Kolmogorov time scale of the turbulence 
$\tau _\eta = (\nu /\varepsilon )^{1/2}$, with 
$\nu$ and 
$\varepsilon$ denoting the kinematic viscosity and the average rate of turbulent dissipation, respectively. Further, the ratios of particle (
$\rho _i$) and fluid (
$\rho _f$) densities are relevant as they characterise to what extent the particle motion is influenced by fluid (‘added mass’) inertia and buoyancy. The particle radii 
$r_i$ determine the collision radius 
$r_c$, and their size relative to the Kolmogorov length scale 
$\eta = (\nu ^3/\varepsilon )^{1/4}$ determines the range of turbulent scales relevant for their motion. In addition to turbulent driving, particle motion may also be affected by gravitational effects, and the relative importance of these two factors is captured by the Froude number 
$Fr = a_\eta /g$, where 
$a_\eta = \eta /\tau _\eta ^2$ and 
$g$ are the turbulence and gravitational accelerations, respectively. Finally, the intensity of the turbulence is measured by the Taylor Reynolds number 
$Re_\lambda = \sqrt {15/(\nu \varepsilon )}\,u'^2$, where 
$u'$ is the single-component root-mean-square (r.m.s.) fluid velocity. Obviously, the entire parameter space spanned by (1.4) is vast and cannot be studied comprehensively here. We therefore limit the present investigation to cases with 
$St_1 = St_2 = St$, and to the zero-gravity regime, i.e. 
$Fr \to \infty$. The benefit of these choices lies in the fact that they keep the problem simple enough to disentangle the relevant mechanisms. Similarly, these configurations are the most amenable to modelling approaches and therefore allow for their evaluation at the most basic level.
Our investigation is based on direct numerical simulations of bubbles and particles in statistically stationary homogeneous isotropic turbulence using the point-particle approach. Details of this approach will be described in § 3 after first reviewing relevant modelling approaches for the collision kernel in § 2. The results are shown in § 4, followed in § 5 by practical considerations in light of our results, and conclusions.
2. Theoretical background and existing models
2.1. The tracer limit 
$St \rightarrow 0$: shear mechanism
In the tracer limit of 
$St \to 0$, the suspended species follow the flow faithfully and distribute uniformly. This means that collisions occur only if particles of finite size are moved relative to each other due to shearing motions in the flow. Considering the dominant shear contribution of the smallest (Kolmogorov) scales of turbulence, and assuming local isotropy as well as Gaussian distribution of the flow velocity gradient, Saffman & Turner (Reference Saffman and Turner1956) derived the classical result
\begin{equation} \varGamma_{12}^{(ST)} = \sqrt{\frac{8{\rm \pi}}{15}}\,\frac{r_c^3}{\tau_\eta}, \end{equation}
predicting the rate of shear-driven collisions in a turbulent flow. Note that here we use the spherical formulation for the collision kernel, which was shown to be the appropriate form by Wang, Ayala & Xue (Reference Wang, Ayala and Xue2005). Employing the concept of a collision cylinder, instead, results in a slightly different value of the prefactor (
$\sqrt {8{\rm \pi} }/3\approx 1.67$ instead of 
$\sqrt {8{\rm \pi} /15} \approx 1.29$).
2.2. Intermediate $St$: preferential sampling and velocity decorrelation
For non-zero 
$St$, the suspended species no longer completely follow the flow. Such inertial effects influence the collision rate via two different pathways. First, even if the drift is small, its accumulated effect leads to preferential concentration, inducing clustering in the particle field. Additionally, the increasing decorrelation between the local particle and fluid velocities affects the collision velocities.
Preferential concentration of inertial particles is widely observed experimentally (Aliseda et al. Reference Aliseda, Cartellier, Hainaux and Lasheras2002; Monchaux, Bourgoin & Cartellier Reference Monchaux, Bourgoin and Cartellier2010; Obligado et al. Reference Obligado, Teitelbaum, Cartellier, Mininni and Bourgoin2014; Petersen, Baker & Coletti Reference Petersen, Baker and Coletti2019; Li et al. Reference Li, Abraham, Guala and Hong2021) and numerically (Bec et al. Reference Bec, Biferale, Cencini, Lanotte, Musacchio and Toschi2007; Goto & Vassilicos Reference Goto and Vassilicos2006; Calzavarini et al. Reference Calzavarini, Kerscher, Lohse and Toschi2008b; Voßkuhle et al. Reference Voβkuhle, Pumir, Lévêque and Wilkinson2014; Ireland, Bragg & Collins Reference Ireland, Bragg and Collins2016). It is rather straightforward to extend (2.1) to account for this effect by simply multiplying it with the RDF 
$g_{12}(r_c)$ (Voßkuhle et al. Reference Voβkuhle, Pumir, Lévêque and Wilkinson2014), yielding
\begin{equation} \varGamma_{12}^{(STc)} = \sqrt{\frac{8{\rm \pi}}{15}}\,\frac{r_c^3}{\tau_\eta}\,g_{12}(r_c). \end{equation}
Various approaches have been proposed to explain the phenomenon of preferential concentration (Maxey Reference Maxey1987; Bec et al. Reference Bec, Celani, Cencini and Musacchio2005, Reference Bec, Biferale, Cencini, Lanotte, Musacchio and Toschi2007; Chen, Goto & Vassilicos Reference Chen, Goto and Vassilicos2006; Goto & Vassilicos Reference Goto and Vassilicos2006; Coleman & Vassilicos Reference Coleman and Vassilicos2009; Zaichik & Alipchenkov Reference Zaichik and Alipchenkov2009; Fouxon Reference Fouxon2012). Among these, the most intuitive one is the ‘centrifuge picture’ (Maxey Reference Maxey1987), according to which heavy particles are ejected out of eddies due to their inertia and hence accumulate in regions of low vorticity and high strain. For collisions between heavy particles (e.g. cloud droplets), it is therefore found that 
$g_{12}(r_c)\geq 1$ generally (Zhou, Wexler & Wang Reference Zhou, Wexler and Wang2001), i.e. clustering enhances the collision rate in some cases even by multiple orders of magnitude (Voßkuhle et al. Reference Voβkuhle, Pumir, Lévêque and Wilkinson2014; Ireland et al. Reference Ireland, Bragg and Collins2016; Pumir & Wilkinson Reference Pumir and Wilkinson2016).
For the collision velocity, inertial effects can be split into a local mechanism and a non-local mechanism. The local mechanism results from the fact that particles react differently to the same fluid forcing provided that they have different properties. Hence this effect contributes an additional relative velocity for bidispersed collisions only, and plays no role in monodispersed cases. For this to occur, the particle trajectories should not deviate significantly from the pathlines (i.e. 
$St$ not too large). An extension of the Saffman–Turner approach accounting for this effect has been reported by Yuu (Reference Yuu1984).
In contrast, the non-local mechanism refers to the situation in which particles arrive at the same location with different particle velocities. Illustratively, one can think of these particles as being ‘slung out’ of neighbouring eddies, and the effect is therefore also known as the ‘sling effect’ (Falkovich, Fouxon & Stepanov Reference Falkovich, Fouxon and Stepanov2002). Unlike the local mechanism, the sling effect is also active for monodisperse suspensions, and so far has been studied almost exclusively in this context (e.g. Falkovich et al. Reference Falkovich, Fouxon and Stepanov2002; Wilkinson, Mehlig & Bezuglyy Reference Wilkinson, Mehlig and Bezuglyy2006; Falkovich & Pumir Reference Falkovich and Pumir2007; Ijzermans, Meneguz & Reeks Reference Ijzermans, Meneguz and Reeks2010; Bewley, Saw & Bodenschatz Reference Bewley, Saw and Bodenschatz2013; Voßkuhle et al. Reference Voβkuhle, Pumir, Lévêque and Wilkinson2014). From these studies, it has become clear (see e.g. Pumir & Wilkinson Reference Pumir and Wilkinson2016, for an overview) that the sling effect can significantly enhance monodisperse collision rates by increasing 
$S_-$. A widely used parametrisation is 
$S_- \sim u_\eta \,F(St,Re_\lambda )$, where 
$u_\eta = \eta /\tau _\eta$ is the Kolmogorov velocity, such that the sling-induced collision rate is given by
\begin{equation} \varGamma_{11}^{(slg)} = 4{\rm \pi} r_c^2 u_\eta\,F(St,Re_\lambda). \end{equation}It has then been proposed (Wilkinson et al. Reference Wilkinson, Mehlig and Bezuglyy2006; Voßkuhle et al. Reference Voβkuhle, Pumir, Lévêque and Wilkinson2014) to obtain the overall collision rate from the sum
\begin{equation} \varGamma_{11}^{(tot)} = \varGamma_{11}^{(STc)} + \varGamma_{11}^{(slg)}. \end{equation}The underlying idea for this decomposition is that particles that are clustered close to each other collide with low velocities, whereas those with high relative velocities can be assumed to be uniformly distributed. Note that both (2.3) and (2.4) are given for the monodisperse case only, since related results for the bidisperse case have not been reported yet.
2.3. Large $St$ limit: kinetic gas behaviour
At very large (but finite) 
$St$, the velocities of the suspended particles arriving at the same point are increasingly uncorrelated. Assuming entirely random and isotropic particle velocities in the spirit of the kinetic gas theory, Abrahamson (Reference Abrahamson1975) derived the collision kernel
\begin{equation} \varGamma_{12}^{(A)} = \sqrt{8{\rm \pi}}\,r_c^2\,\sqrt{v_1'^2 + v_2'^2}, \end{equation}
where the mean-square single-component particle velocity 
$v_i'^2$ is related to flow properties via
\begin{equation} v_i'^2 = A_iu'^2 = \frac{\dfrac{T_{fL}}{\tau_i}+\gamma_{i}^2}{\dfrac{T_{fL}}{\tau_i}+1}\,u'^2, \end{equation}
with 
$T_{fL}$ denoting the fluid Lagrangian integral time scale, and 
$\gamma _i=3\rho _f/(2\rho _i + \rho _f)$. It has been pointed out (Voßkuhle et al. Reference Voβkuhle, Pumir, Lévêque and Wilkinson2014) that (2.5) is not strictly valid for turbulence as it fails to account for the multiscale structure of the flow. Alternative derivations (Völk et al. Reference Völk, Jones, Morfill and Röser1980; Mehlig, Uski & Wilkinson Reference Mehlig, Uski and Wilkinson2007) based on the Kolmogorov (Reference Kolmogorov1941) phenomenology arrived at 
$F(St,\infty ) \sim K\sqrt {St}$, with 
$K$ being a universal dimensionless constant, at the limit of intense turbulence in the context of (2.3).
2.4. Modelling approaches for bubble–particle collisions in turbulence
Next, we outline briefly the different approaches to modelling bubble–particle collisions reported in particular in the mining literature so far. For additional details, we refer the reader to the reviews on the topic (e.g. Nguyen et al. Reference Nguyen, An-Vo, Tran-Cong and Evans2016; Hassanzadeh et al. Reference Hassanzadeh, Firouzi, Albijanic and Celik2018; Kostoglou, Karapantsios & Oikonomidou Reference Kostoglou, Karapantsios and Oikonomidou2020b). Note that here and in the following we use subscripts 
$i = b,p$ to denote bubbles and heavy particles, respectively.
The most commonly adopted approach is to assume the high-
$St$ limit and to base the collision models on (2.5). The differences from the theory of Abrahamson (Reference Abrahamson1975) are related to the expressions used to model the r.m.s. velocities. Instead of (2.6), Schubert (Reference Schubert1999) and later Bloom & Heindel (Reference Bloom and Heindel2002) used the relation given by Liepe & Möckel (Reference Liepe and Möckel1976), which reads
\begin{equation} w_i'^{(LM)} = 0.57\,\frac{\varepsilon^{4/9}r_i^{7/9}}{\nu^{1/3}}\left(\frac{\rho_f-\rho_i}{\rho_f}\right)^{2/3}. \end{equation}
This result is based on an analogy to gravitational settling with the fluid acceleration in the inertial range replacing gravity. The resulting apparent weight is balanced by Allen's drag, which scales with the particle slip velocity 
$\boldsymbol {w}_i = \boldsymbol {v}_i - \boldsymbol {u}$ as 
$|\boldsymbol {w}_i|^{3/2}$ (with boldface denoting vectors). Here, 
$\boldsymbol {v}_i$ is the bubble/particle velocity, and 
$\boldsymbol {u}$ is the flow velocity. Later work by Nguyen & Schulze (Reference Nguyen and Schulze2004) used (2.7) with a different constant, 0.83, for bubbles. For particles, they replaced the inertial subrange acceleration with that in the dissipation range, and Allen's drag with Stokes drag, in order to account for the small size of typical particles. This resulted in
\begin{equation} w_p'^{(NS)} = \frac{2r_p^3\varepsilon}{135\nu^2}\left(\frac{\rho_p-\rho_f}{\rho_f}\right). \end{equation}
Importantly, both (2.7) and (2.8) are expressions for the relative velocity between bubbles/particles and the surrounding fluid. Their use with Abrahamson's theory is therefore inconsistent since (2.5) contains velocities in a fixed frame of reference. This was already pointed out in Kostoglou et al. (Reference Kostoglou, Karapantsios and Evgenidis2020a). We further note that the quasi-static assumption underlying (2.7) and (2.8) is equivalent to a low-
$St$ approximation and therefore not valid for the high-
$St$ limit in which (2.5) applies. In fact, (2.8) is consistent with the rigorously derived small-
$St$ limit (Fouxon Reference Fouxon2012)
\begin{equation} \boldsymbol{w}_i=\beta_i\tau_i \boldsymbol{a}_f, \end{equation}with
\begin{equation} \beta_i = \frac{2(\rho_f - \rho_i)}{\rho_f + 2\rho_i}, \end{equation}
if the fluid acceleration is approximated by the dissipative scaling 
$|\boldsymbol {a}_f| \approx \varepsilon r_i/(15 \nu )$. Note, however, that it appears more appropriate to use 
$|\boldsymbol {a}_f| \approx a_\eta$ for small particles (
$r_i/\eta \lessapprox 1$) because 
$|\boldsymbol {a}_f| \rightarrow 0$ otherwise. In the general case, the fact that the fluid acceleration experienced by the particle may vary considerably over the particle response time precludes a simple relation between 
$\boldsymbol {w}_i$ and 
$\boldsymbol {a}_f$.
In a different approach, Ngo-Cong, Nguyen & Tran-Cong (Reference Ngo-Cong, Nguyen and Tran-Cong2018) extended the model by Yuu (Reference Yuu1984) to the bubble–particle case. The resulting expression takes the form
\begin{equation} \varGamma_{bp}^{(NC)} = (r_b + r_p)^2\sqrt{\frac{8{\rm \pi}}{3}\left\{3(A_b + A_p - 2B)u'^2 + (A_br_b^2 + A_pr_p^2 + 2Br_pr_b)\,\frac{\varepsilon}{3\nu}\right\}}, \end{equation}
where the term proportional to 
$u'^2$ represents the ‘local’ inertial effect on the relative velocity that adds to the shear-driven collisions that are accounted for by the 
$\varepsilon$ term. Aside from bubble/particle properties, the coefficients 
$A_i$ and 
$B$ depend also on 
$T_{fL}$, similar to (2.6). The bubble–particle velocity correlation is determined by 
$B$. Ngo-Cong et al. (Reference Ngo-Cong, Nguyen and Tran-Cong2018) additionally incorporated (2.7) and (2.8) into this model. However, doing so suffers from the same inconsistencies outlined above. As a consequence, this expression featured a negative radicand when evaluated for the parameters in this study, and is therefore not included further. It is worth mentioning that the term proportional to 
$u'^2$ in (2.11) does not approach (2.5) even at large 
$St$, since 
$A_i=B$ for identical particles. We have therefore extended the model along the lines of the theory of Kruis & Kusters (Reference Kruis and Kusters1997), which attempted to reconcile this issue by marrying the concept of Yuu (Reference Yuu1984) for the small 
$St$ case and that of Williams & Crane (Reference Williams and Crane1983) (which does not consider shear-driven collisions) for the large 
$St$ case, to bubble–particle collisions (see Appendix A). The original formulation of Kruis & Kusters (Reference Kruis and Kusters1997) accounted for collisions of particles of equal density only. We further note the work of Fayed & Ragab (Reference Fayed and Ragab2013), who employed the model of Zaichik, Simonin & Alipchenkov (Reference Zaichik, Simonin and Alipchenkov2010), which was developed originally for collisions between particles with arbitrary but equal density assuming a joint-normal fluid–particle velocity distribution, to the bubble–particle case. Recently, Kostoglou et al. (Reference Kostoglou, Karapantsios and Evgenidis2020a) proposed another model that is specific to the case of very fine particles that essentially follow the flow, such that the relative velocity is dominated by the bubble slip velocity.
A common feature of almost all these models is that they are direct adaptations of concepts developed for collisions between heavy particles. They therefore fail to account for fundamental differences in how bubbles and heavy particles react to a turbulent flow. Most strikingly, for example, the response to fluid accelerations given in (2.9) is in opposite directions as 
$\beta$ switches sign from negative to positive from 
$\rho _i>\rho _f$ to 
$\rho _i<\rho _f$. Similarly, applying the centrifuge picture to light particles such as bubbles, one expects them to concentrate in regions of high vorticity as they travel to the centre of the eddies. This implies that heavy and light particles segregate in a turbulent flow, as has indeed been observed (Calzavarini et al. Reference Calzavarini, Cencini, Lohse and Toschi2008a; Fayed & Ragab Reference Fayed and Ragab2013). As a consequence, it is expected that 
$g_{bp}(r_c)< 1$ in these cases, such that preferential concentration is expected to lead to a decrease in 
$\varGamma _{bp}^{(STc)}$ for collisions between heavy and light particles. Such aspects cannot be quantified easily from laboratory-scale flotation set-ups (e.g. Darabi et al. Reference Darabi, Koleini, Deglon, Rezai and Abdollahy2019) since it is difficult to disentangle the many factors influencing the overall flotation rate. Numerically, Wan et al. (Reference Wan, Yi, Wang, Sun, Chen and Wang2020) observed no reduction of the RDF, but increased relative velocities for bubble–particle pairs based on a point-particle simulation while taking gravitational effects into account. For cases without gravity, Fayed & Ragab (Reference Fayed and Ragab2013) reported 
$g_{bp}(r_c) < 1$ and relative velocities matching those predicted by the theory of Zaichik et al. (Reference Zaichik, Simonin and Alipchenkov2010) for simulations of the extreme case 
$\rho _b/\rho _f = 0$ and 
$\rho _p/\rho _f\to \infty$ across a range of 
$\tau _i$. It is the goal of the present study to add to this a systematic investigation of how bubble–particle relative behaviour affects their collision statistics in the low- to moderate-
$St$ regime, and to assess how this affects modelling outcomes.
3. Methods
3.1. Fluid phase
To obtain the background turbulence, we solve the Navier–Stokes equation and the continuity equation for incompressible flow:
$$\begin{gather} \frac{\mathrm{D}\boldsymbol{u}}{\mathrm{D}t} ={-}\frac{1}{\rho_f}\,\boldsymbol{\nabla} P + \nu\, \nabla^2\boldsymbol{u} + \boldsymbol{f}_\varPsi, \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0, \end{gather}$$
where 
$\mathrm {D}/\mathrm {D}t$ is the material derivative following a fluid element, 
$t$ is the time, and 
$P$ is the pressure. The forcing 
$\boldsymbol {f}_\varPsi$, which is non-zero only for the wavenumbers 
$|\boldsymbol {\kappa }|/|\boldsymbol {\kappa }_0| < 2.3$ (i.e. the largest scales), with 
$|\boldsymbol {\kappa }_0|$ being the smallest wavenumber along each direction, is added to counter dissipation and maintain statistical stationarity. We employ the widely used Eswaran & Pope (Reference Eswaran and Pope1988) forcing scheme (Chouippe & Uhlmann Reference Chouippe and Uhlmann2015; Spandan et al. Reference Spandan, Putt, Ostilla-Mónico and Lee2020). In brief, a complex vector is generated in Fourier space for the forced wavenumbers by multiple Uhlenbeck–Ornstein processes. This vector is then projected onto the plane normal to the wavevector, thereby ensuring that 
$\boldsymbol {f}_\varPsi$ is divergence-free. To simulate the fluid motion, a second-order finite-difference solver is implemented on a staggered grid (Verzicco & Orlandi Reference Verzicco and Orlandi1996; van der Poel et al. Reference van der Poel, Ostilla-Mónico, Donners and Verzicco2015). All spatial derivatives, including the nonlinear terms, are discretised by a second-order central finite-difference method. The simulation domain is a cubic box with length 
$L_{box}=1$ and triply periodic boundary conditions to eliminate boundary effects. Time marching is performed using a fractional step third-order low-storage Runge–Kutta scheme and the implicit Crank–Nicolson scheme for all viscous terms at a maximum Courant–Friedrichs–Lewy (CFL) number 1.2. The CFL number is 
$\max (|u_1| + |u_2| + |u_3|)\,\Delta t/\Delta x$ over all cells, where the grid spacing 
$\Delta x$ is identical in all the three dimensions. Here, 
$u_1,u_2,u_3$ are the 
$x$-, 
$y$-, 
$z$-components of 
$\boldsymbol {u}$, and 
$\Delta t$ is the time step. The simulation is parallelised via slab decomposition along the 
$z$-direction.
Turbulence with 
$Re_\lambda = 72$ and 
$175$ is generated from an initially quiescent fluid following the tuning method proposed by Chouippe & Uhlmann (Reference Chouippe and Uhlmann2015). The simulations are allowed to run until the pseudo-dissipation 
$\bar {\varepsilon }$ and 
$Re_\lambda$ are statistically stationary over at least 
$30T_L$. These statistically stationary flow fields are then used as the starting fields for the point-particle simulations. The flow statistics are listed in table 1. For validation, the longitudinal and transverse energy spectra are plotted in figure 1. Excellent agreement with the literature is shown.
Statistics of the homogeneous isotropic turbulence: the grid size (
$\mathcal {N}$), pseudo-dissipation (
$\bar {\varepsilon }$), Kolmogorov length scale (
$\eta$), maximum wavenumber (
$k_{max}$), Kolmogorov velocity (
$u_\eta$) scale, r.m.s. velocity fluctuations (
$u'$), large-scale isotropy (
$u_x'/u_y'$), and large eddy turnover time (
$T_L$) relative to the Kolmogorov time scale (
$\tau _\eta$). Here, 
$N_b,N_p$ are the numbers of bubbles and particles, respectively.
The (a) longitudinal and (b) transverse energy spectra in single-phase statistically stationary homogeneous isotropic turbulence. Dashed lines show the data from Jiménez et al. (Reference Jiménez, Wray, Saffman and Rogallo1993); triangles represent the 
$-5/3$ power law. The agreement with the literature is excellent, so the dashed lines can be obscured, especially in (b).
3.2. Suspended phases
Bubbles and particles in the system are modelled using the point-particle approximation, where forces act on point masses. The bubble and particle dynamics are governed by (Maxey & Riley Reference Maxey and Riley1983; Tchen Reference Tchen1947)
\begin{equation} \frac{4}{3}\,{\rm \pi} r_i^3\rho_i\,\frac{{\rm d}\boldsymbol{v}_i}{{\rm d}t} = 6{\rm \pi}\mu r_i\,f_i(\boldsymbol{u}-\boldsymbol{v}_i) + \frac{4}{3}\,{\rm \pi} r_i^3\rho_f\,\frac{\mathrm{D}\boldsymbol{u}}{\mathrm{D}t} + \frac{2}{3}\,{\rm \pi} r_i^3\rho_f\left(\frac{\mathrm{D}\boldsymbol{u}}{\mathrm{D}t} - \frac{{\rm d}\boldsymbol{v}_i}{{\rm d}t}\right), \end{equation}
where 
$\mu =\nu \rho _f$ is the absolute viscosity, 
$r_i$ is determined from 
$St_i$, and 
$\boldsymbol {u}$ is evaluated at particle position in this equation. The three terms on the right-hand side of (3.2) are the drag force, the pressure gradient force and the added mass force, respectively. Note that history forces, lift and reverse coupling are not considered to render it easier to disentangle the relative behaviour of bubble and particles. Nevertheless, as the lift force can be expected to play a role for bubbles, its effect will be discussed in § 4.3. Furthermore, a one-way coupled system is a consistent choice to study the geometric collision rate where hydrodynamic interactions are neglected. It entails the assumption of the dilute limit in which turbulence modifications by the suspended species are negligible (Brandt & Coletti Reference Brandt and Coletti2022). The correction factor 
$f_i = 1 + 0.169\,Re_i^{2/3}$ accounts for finite bubble or particle Reynolds number 
$Re_i = 2r_i\,|\boldsymbol {w}_i|/\nu$, and implies the assumption of rigid spheres that obey the no-slip boundary condition for both species (Nguyen & Schulze Reference Nguyen and Schulze2004). This is realistic since liquids in flotation cells typically contain significant amounts of surfactants, so that the bubble surfaces would likely be contaminated (Nguyen & Schulze Reference Nguyen and Schulze2004; Huang, Legendre & Guiraud Reference Huang, Legendre and Guiraud2012). Although other commonly used expressions for 
$f_i$ are available, the difference is minimal, as shown in Appendix C.
Equation (3.2) is solved for each bubble and particle using a finite-difference scheme. To determine the flow velocities and the velocity gradients at bubble and particle positions, these quantities are interpolated from the Eulerian grid of the fluid solver to the particle positions using tri-cubic Hermite spline interpolation with a stencil of four points per direction. This choice is made as Hermite splines are comparable in accuracy (van Hinsberg, Clercx & Toschi Reference van Hinsberg, Clercx and Toschi2017) and computationally cheaper to implement than B-splines (Ostilla-Monico et al. Reference Ostilla-Monico, Yang, van der Poel, Lohse and Verzicco2015). Time marching of 
$\boldsymbol {v}_i$ is performed with the explicit forward Euler method, and that of the positions of the suspended phases is done using the second-order Adams–Bashforth scheme. For stability, the time step is restricted such that neither the fluid nor the particle CFL number (
$\max (v_{i1},v_{i2},v_{i3})\,\Delta t/\Delta x$) exceeds the value 1.2. This limit is enforced in both the suspended- and fluid-phase solvers, which run with the same simulation time step. We have compared particle statistics to data from the literature in order to verify our code, and those results are included in Appendix B.
Collisions between particles and bubbles are treated as ‘ghost collisions’. Under this approach, a ‘collision’ occurs once the centres of the members of an approaching pair reach the collision distance and the colliding pair pass each other without interaction. The collision radii are determined through the virtual radii 
$r_i$ as computed from 
$St_i$. This scheme is often employed by simulations of particles in turbulence (Bec et al. Reference Bec, Celani, Cencini and Musacchio2005, Reference Bec, Biferale, Cencini, Lanotte, Musacchio and Toschi2007; Goto & Vassilicos Reference Goto and Vassilicos2008; Voßkuhle et al. Reference Voβkuhle, Pumir, Lévêque and Wilkinson2014; Ireland et al. Reference Ireland, Bragg and Collins2016) and has been shown to be consistent with the formulation of 
$\varGamma _{pp}^{(ST)}$ (Wang, Wexler & Zhou Reference Wang, Wexler and Zhou1998). To suppress the effect of different length scales when comparing collision pairs, we take 
$r_c = r_b + r_p$ for every type of collision. This is unlike studies examining solely particle–particle collisions, which define 
$r_c = 2r_p$ (e.g. Voßkuhle et al. Reference Voβkuhle, Pumir, Lévêque and Wilkinson2014; Ireland et al. Reference Ireland, Bragg and Collins2016). Numerically, these collisions are detected with the ‘proactive’ detection scheme in Sundaram & Collins (Reference Sundaram and Collins1996).
The details of the simulation of the suspended phases are as follows. Particles (
$\rho _p/ \rho _f = 5$) and bubbles (
$\rho _b/ \rho _f = 1/1000$) with 
$St_i = 0.1\unicode{x2013}3$ were seeded randomly and homogeneously (
$10\,000\unicode{x2013}140\,000$ each) into the turbulent flow after it had reached a statistically stationary state as described in § 3.1. The density ratios correspond to sulphide minerals colliding with air bubbles in water, while the simulated 
$St$ range falls within 
$0.1 \lesssim St \lesssim 100$, which corresponds to particle sizes yielding the highest mineral recovery rate in conventional flotation cells (
$10\,\mathrm {\mu }{\rm m}\lesssim r_p\lesssim 100\,\mathrm {\mu }{\rm m}$ when 
$\nu = 10^{-6}\,{\rm m}^{2}\,{\rm s}^{-1}$ and 
$\varepsilon = 5\,{\rm m}^{2}\,{\rm s}^{-3}$, following Ngo-Cong et al. Reference Ngo-Cong, Nguyen and Tran-Cong2018). For each case, bubbles and particles have identical 
$St_i$ to focus on the effect of different densities and to keep the parameter space manageable. The upper limit 
$St_i \leq 3$ is mandated by the fact that the increasingly large (virtual) bubble radius for larger 
$St$ violates the point-particle approximation. For 
$St_b = 3$ and at 
$Re_\lambda = 175$, we have 
$r_b/\eta \approx 5$, which is already marginal (Homann & Bec Reference Homann and Bec2010). We monitor the p.d.f.s of bubble and particle positions for statistical stationarity. Once this has been reached, collision statistics are collected over at least 
$7.7T_L$. Two kinds of statistics are evaluated: in a fixed reference frame considering all particles/bubbles (Eulerian statistics), and along individual trajectories of the pairs that collide (Lagrangian statistics). The former are computed at least every 
$\sim 0.06T_L$, while the latter are calculated every time step (with only 25 % of all bubbles for the bubble–bubble case as the number of colliding bubbles is high). All the statistics presented are time- and ensemble-averaged unless indicated otherwise.
4. Results
4.1. Eulerian statistics
4.1.1. Collision kernel
Figure 2(a) shows simulation results for the dimensionless bubble–particle (bp), bubble–bubble (bb) and particle–particle (pp) collision kernels 
$\varGamma$. In addition to determining 
$\varGamma$ based on counting the number of collisions per time step (solid triangles), the collision kernels for 
$Re_\lambda = 175$ are determined indirectly via the RDF and the effective radial collision velocity according to (1.2) (shown as dots). Both results match closely, verifying our analysis procedure. With the normalisation by 
$\tau _\eta /r_c^3$ suggested by the Saffman–Turner framework, the results are insensitive to the change in 
$Re_\lambda$ for the pp and bp cases, but not for bb collisions. From the data, it is further evident that the relative behaviour of bubbles and particles is distinct from that of identical particles. For collisions between identical species, 
$\varGamma \tau _\eta /r_{c}^3$ is maximum when 
$St \sim 1$, while 
$\varGamma _{bp}\tau _\eta /r_{c}^3$ exhibits a minimum for this value of 
$St$. This trend is not captured by any of the models for the bp case discussed in § 2.4, for which the predictions are included as lines in figure 2(a). Generally, the model predictions are also significantly higher than the actual collision rates obtained from the simulations, the exception being the Saffman & Turner (Reference Saffman and Turner1956) model, which best captures the magnitude yet fails to predict the proper 
$St$-trend, as is illustrated more clearly in figure 2(b).
(a) The dimensionless bubble–particle (bp), bubble–bubble (bb) and particle–particle (pp) collision kernels at collision distance. The triangle symbols denote 
$\varGamma$ determined directly from the collision rate, while 
$\varGamma _{4{\rm \pi} }$ (shown only for the 
$Re_\lambda = 175$ cases) is calculated according to (1.2). Unless specified otherwise, the colour conventions for all figures follow this figure. (b) Zoomed-in version of (a) plotted on linear scale.
The fact that the Abrahamson (Reference Abrahamson1975) and large-
$St$ Kruis & Kusters (Reference Kruis and Kusters1997) predictions do not match the data is to be expected as the present 
$St$ range does not match the assumptions made in these frameworks. Naturally, this also transfers to all approaches based on 
$\varGamma _{bp}^{(A)}$, and the somewhat better agreement with our data for models employing (2.7) and (2.8) instead of (2.6) is rather an artefact of the inconsistencies discussed earlier. This is emphasised by figure 3(a), where the large difference between the modelled slip velocities and the mean-square bubble/particle velocities 
$v_i'^2$ is obvious. Notably from the same figure, 
$v_p'^2$ is well predicted by the models of Abrahamson (Reference Abrahamson1975), Kruis & Kusters (Reference Kruis and Kusters1997) (small 
$St$) and Zaichik et al. (Reference Zaichik, Simonin and Alipchenkov2010). On the other hand, model estimates for 
$v_b'^2$ are generally too high. It is insightful to note the stark overprediction by (2.6) as the underlying framework by Abrahamson (Reference Abrahamson1975) is also employed in the models by Ngo-Cong et al. (Reference Ngo-Cong, Nguyen and Tran-Cong2018). The resulting overprediction of 
$v_b'^2$ might therefore explain in part why the estimate of 
$\varGamma _{bp}^{(NC)}$ is too large even though their concept (which is based on Yuu (Reference Yuu1984)) is in principle more suitable for the moderate values of 
$St$ here. Unlike Abrahamson (Reference Abrahamson1975), Kruis & Kusters (Reference Kruis and Kusters1997) (small 
$St$) include dissipation range scaling in their modelling, which is seen to improve the prediction at 
$St\lessapprox 1$ but still overestimates 
$v_b'^2$ for larger 
$St$. Finally, there is another factor, which affects all models. Bubbles especially do not sample the flow randomly, such that the mean-square fluid velocity at bubble locations 
$\langle u'^2\rangle _b$ is almost 10 % lower than 
$u'^2$, as shown in figure 3(b). This implies that bubbles sample flow regions with weaker fluid velocity fluctuations such that using 
$u'^2$ as model input may overestimate 
$v'^2_b$. Figure 3(b) also shows that this effect is less pronounced for heavy particles. The underlying cause for these observations is preferential sampling of flow regions, which we will discuss in § 4.1.2.
(a) The mean-square velocity of bubbles and particles at various 
$St$. Models are shown for 
$Re_\lambda = 175$ only. (b) The mean-square fluid velocity conditioned at bubble/particle positions for 
$Re_\lambda = 175$.
4.1.2. Bubble/particle spatial distribution
In order to elucidate in particular the 
$St$ trends for 
$\varGamma _{bp}$, we first investigate the distribution of bubbles and particles in the flow. While all the models assume homogeneous distributions, this does not hold at intermediate 
$St$, as figure 4 confirms, where instantaneous snapshots of the bubble and particle fields at different 
$St$ and 
$Re_\lambda =175$ are shown. For 
$St = 1$ and 3, bubbles and particles are seen to cluster but do so in different regions of the flow. This behaviour has been observed previously in the literature (Calzavarini et al. Reference Calzavarini, Cencini, Lohse and Toschi2008a; Fayed & Ragab Reference Fayed and Ragab2013; Wan et al. Reference Wan, Yi, Wang, Sun, Chen and Wang2020) and is additionally shown by the different mean-square fluid velocity at bubble/particle positions in figure 3(b). To investigate this preferential concentration, we plot the norm of the rotation 
$\langle R^2 \rangle _r$ and strain 
$\langle \mathcal {S}^2 \rangle _r$ rates of the flow at the particle/bubble positions in figure 5, with 
$\langle \cdot \rangle _r$ denoting an ensemble average over particles in pairs with separations smaller than 
$r$. For tracers (
$St=0$), we obtain 
$\tau _\eta ^2\langle R^2 \rangle _r=\tau _\eta ^2\langle \mathcal {S}^2 \rangle _r =0.5$ with 
$r\to \infty$, consistent with the analytical result in statistically stationary homogeneous isotropic turbulence. For the other cases, 
$\langle R^2\rangle _r$ and 
$\langle \mathcal {S}^2\rangle _r$ are conditioned on pairs with 
$r\leq 2r_c$. Figure 5(a) shows that bubbles (particles) cluster in regions of high (low) rotation rate, which is consistent with the centrifuge picture. Due to the clustering, conditioning has little effect for monodisperse collisions, and the results are therefore very close to single particle statistics. This is different for bp pairs, where 
$R^2$ is consistently lower (higher) for bubbles (particles) close to a particle (bubble). These observations imply that bp collisions occur for a subset of bubbles/particles that is located outside of their respective ‘preferred’ location within the flow. We note that correlations between particle/bubble locations and the strain rate (figure 5b) are much weaker and do not display the same qualitative trends observed for 
$R^2$. This is consistent with observations by e.g. Ireland et al. (Reference Ireland, Bragg and Collins2016) and Wang et al. (Reference Wang, Wan, Yang, Wang and Chen2020), who also found only a weak correlation between the positions of heavy particles and regions of the flow with low strain rate.
Instantaneous snapshots of bubbles and particles in a slice with width 
$\times$ height 
$\times$ depth = 
$L_{box} \times L_{box} \times 20\eta$ in the non-transient state for 
$Re_\lambda = 175$ at (a) 
$St = 0.1$, (b) 
$St = 1$, and (c) 
$St = 3$. The size of the 
$St = 0.1$ bubbles and particles is tripled for visibility.
Average value of the norm of (a) the rotation rate 
$R^2$, and (b) the strain rate 
$\mathcal {S}^2$, of the flow at bubble/particle/tracer positions. The bubble and particle data are conditioned on pairs with separation 
$r \leq 2r_{c}$. Here, bp,b (bp,p) refer to bubbles (particles) in bubble–particle pairs.
A consequence of the preferential concentration in different flow regions is that bubbles and particles become segregated. This effect is quantified by the RDF 
$g(r)$. Essentially, 
$g(r)$ relates the actual number of pairs with separation 
$r$ to that expected for uniformly distributed particles. Therefore, 
$g(r) = 1$ when particles are uniformly distributed, 
$g(r)>1$ implies clustering, and 
$g(r)<1$ implies segregation (Saw Reference Saw2008). Figure 6(a) shows results for the RDF at collision distance from our simulations. Bubble–particle pairs indeed segregate, and the segregation is strongest when 
$St = 1$, as reflected by the local minimum of 
$g(r_{c})$. This is compatible with the corresponding behaviour of bubbles and particles, which always cluster for the tested parameters, with maximum clustering when 
$St\sim 1$. Increasing 
$Re_\lambda$ increases segregation slightly, but the effect is weak.
(a) RDF at collision distance. (b) The bubble–particle RDF as a function of 
$r$ and the segregation length scale 
$r_{seg}$ at 
$Re_\lambda = 175$. The number above each line is the corresponding 
$St$, and the dotted segments denote the RDF below the collision distance. Inset shows 
$r_{seg}$ at various 
$St$.
Figure 6(b) displays the bubble–particle RDF as a function of 
$r$ at 
$Re_\lambda =175$, which helps us to understand the effect of different length scales. The shape of the 
$g(r)$ curves is distinctly different from the power-law behaviour reported for the monodisperse case (Ireland et al. Reference Ireland, Bragg and Collins2016). For all 
$St$, there is an essentially flat region for 
$g(r)$ at small scales, followed by a transition region with the steepest gradient at intermediate scale 
${\sim }10\eta$ before approaching 1 for large separations. This behaviour suggests that the relevant length scale for the segregation, i.e. a typical distance between bubble and particle clusters, is at intermediate scales. To quantify this more precisely, we define the separation corresponding to the point of inflexion in 
$g_{bp}(r)$ as the segregation length scale 
$r_{seg}$. The values of 
$r_{seg}$ are marked by crosses in figure 6(b), and plotted against 
$St$ in the inset. The figure shows that 
$r_{seg}/\eta$ increases approximately linearly with 
$St$ and increases slightly with 
$Re_\lambda$. These findings are consistent with those by Calzavarini et al. (Reference Calzavarini, Cencini, Lohse and Toschi2008a), who studied segregation based on a concept inspired by Kolmogorov's distance measure (Kolmogorov Reference Kolmogorov1963). In particular, these authors also report segregation scales 
${\sim }10\eta$ with an increasing trend for higher 
$Re_\lambda$.
The trends observed for 
$g(r_c)$ in figure 6(a) remarkably resemble those discussed for the 
$St$ dependence of 
$\varGamma$ earlier, in the context of figure 2. It is therefore suggestive to think that the effects of segregation may explain in particular the discrepancy between the data and 
$\varGamma ^{(ST)}$. We can check this by plotting 
$\varGamma ^{(STc)}$ as defined in (2.2), which is shown as hollow symbols in figure 7(a). From this plot, it can be seen that correcting 
$\varGamma ^{(ST)}$ with the RDF leads to an almost perfect match with the pp collision kernel. This differs from, but is only seemingly at odds with, the results in Voßkuhle et al. (Reference Voβkuhle, Pumir, Lévêque and Wilkinson2014) due to the larger collision distance considered here as a consequence of keeping 
$r_c= r_b + r_p$ constant for all collision types. This therefore implies that while pp collisions at moderate 
$St$ are governed by the sling mechanism at small separations, they remain shear-dominated at larger ones. In contrast, 
$\varGamma ^{(STc)}$ overcorrects the bp collision kernel and undercorrects that of bb collisions. In figure 7(b), we replot the same data in the form of the ratio 
$\varGamma ^{(STc)}/\varGamma$, which can be interpreted as the relative contribution of the shear mechanism (compensated for segregation/clustering) to the overall collision rate. For pp collisions, this ratio is very close to 1 throughout. The situation is different for the bb case, where 
$\varGamma ^{(STc)}/\varGamma$ approaches 1 only for the lowest 
$St$ considered, and the value drops significantly for the higher 
$St$. Consistently the lowest values for 
$\varGamma ^{(STc)}/\varGamma$ are observed for bp collisions where the value quickly drops to around 0.5, which implies that a significant part of the relative velocities cannot be explained by the shear mechanism in this case. It is important to stress here that keeping 
$r_c$ constant for all collision types for a given 
$St$ removes this as a factor, such that the observed trends across species can only be rooted in differences in the relative approach velocities, which will be studied next.
(a) The dimensionless collision kernel (solid symbols) and the Saffman–Turner prediction after accounting for the RDF 
$\varGamma ^{(STc)}$ (hollow symbols). (b) The ratio between 
$\varGamma ^{(STc)}$ and the collision kernel.
4.1.3. Relative velocities
We show results for the non-dimensional effective radial approach velocity 
$S_-/u_\eta$ at collision distance in figure 8(a). The non-dimensional approach velocity increases with 
$St$ and also slightly with 
$Re_\lambda$ for collisions involving bubbles, whereas it has only a very weak 
$Re_\lambda$ dependence for pp collisions, as also reported in Ireland et al. (Reference Ireland, Bragg and Collins2016). Consistent with figure 7, the approach velocities are highest for bp collisions and lowest in the pp case.
(a) The effective radial collision velocity. (b) The effective bubble–particle radial approach velocity against pair separation 
$r$ at 
$Re_\lambda = 175$. (c) The effective radial approach velocity of each type of collisions across 
$St = 0.1$ to 
$St = 3$ at 
$Re_\lambda = 175$. The dashed lines show the prediction by Zaichik et al. (Reference Zaichik, Simonin and Alipchenkov2010).
To understand these trends, it is instructive to consider 
$S_-$ as a function of 
$r$, which is presented in figure 8(b) for bp collisions for different 
$St$, and in the plots of figure 8(c) for different collision types at constant 
$St$. Also shown in these figures is the effective approach velocity for tracer particles, for which 
$S_-$ is proportional to 
$r$ in the dissipation range, as expected. In particular, for low 
$St$, where 
$S_-(r)$ remains close to the tracer curve, this in part explains the 
$St$ dependence of 
$S_-(r_c)$ as 
$r_c$ increases with increasing 
$St$. However, the curves also ‘peel off’ the linear scaling at increasingly larger distances and to a larger extent as 
$St$ increases. For the pp case, this is a well-known manifestation of the ‘sling effect’, where due to path history effects, particles arrive at the same location with different velocities such that their relative velocities exceed those of the fluid. Note that while the bb case has a larger 
$S_-$ than the tracer case, the linear scaling remains mostly unchanged there, such that the difference presumably is rather due to preferential concentration effects and the fact that the slip velocity at the same 
$St$ is higher for this case (see (2.9)). The deviations from the tracer behaviour are strongest for the bp case throughout, and already for 
$St\geq 2$ 
$S_-$ are essentially independent of 
$r$ for 
$r\lessapprox 10\eta$. The decorrelation between bubble and particle velocities implicit in these observations is the reason why 
$S_-$ is largest for all 
$St$ in figure 8(a). The model of Zaichik et al. (Reference Zaichik, Simonin and Alipchenkov2010) fails to reproduce this feature. Instead, 
$S_-$ from the model is highest for bb collisions, such that the trend falls in line with that of 
$v_i'^2$ (see figure 3). Plotting the model prediction against the separation distance in figure 8(c) confirms that this trend persists over all 
$r$. Approach velocities are also generally overpredicted by the model, which is in part also related to the differences observed for 
$v_b'^2$ in figure 3 and to the significant peel off the linear scaling especially for the bb case. Consistent with the results in figure 7, 
$S_-^{pp}$ at both 
$Re_\lambda$ is matched very well by the Saffman–Turner model 
$S_-^{(ST)} = 0.5r_c\sqrt {2\varepsilon /15\nu {\rm \pi}}$, where the prefactor 0.5 accounts for the fact that only separating pairs are considered here.
In modelling approaches (e.g. Abrahamson Reference Abrahamson1975; Yuu Reference Yuu1984; Kruis & Kusters Reference Kruis and Kusters1997; Zaichik et al. Reference Zaichik, Simonin and Alipchenkov2010; Ngo-Cong et al. Reference Ngo-Cong, Nguyen and Tran-Cong2018), it is common to deduce 
$S_-$ based on a Gaussian distribution of the relative radial velocity 
$\Delta v_r$. In this case, the ratio 
$S_-/\sqrt {S_{2\parallel }}$ can be determined to be 
$1/\sqrt {2{\rm \pi} } \approx 0.4$, where the variance is defined as
\begin{equation} S_{2\parallel}(r) = \int_{-\infty}^{\infty} (\Delta v_r)^2\mathrm{p.d.f.}(\Delta v_r|r)\,\mathrm{d}(\Delta v_r). \end{equation}
It was already reported in Ireland et al. (Reference Ireland, Bragg and Collins2016) that 
$S_-/\sqrt {S_{2\parallel }}$ can drop significantly below the Gaussian value for pp collisions at 
$St \sim 1$ and for small 
$r_c$. This is confirmed by our results in figure 9, where also the bb case is seen to follow similar trends. For bp collisions, the ratio remains much closer to the Gaussian value, almost indifferent to 
$St$ and close to the tracer result. Relative velocities are therefore better approximated by the Gaussian assumption for the bp case, while the overprediction of 
$S_-$ resulting from doing so is larger for the monodisperse collisions.
The ratio of the effective radial approach velocity to the standard deviation of the radial component of the relative velocity 
$\Delta v_r$ at 
$Re_\lambda = 175$.
For a more quantitative investigation of variation of 
$S_-$ across collision types, we plot the excess bubble–particle effective radial collision velocity 
$\Delta S_- = S_-^{bp} - (S_-^{bb}+S_-^{pp})/2$ in figure 10(a). At small 
$St$, such an additional approach velocity arises for the bp case from the fact that, due to the change in sign of 
$\beta _i$ in (2.9), bubbles and particles react differently when subjected to the same fluid acceleration. Indeed, figure 10(b) shows that 
$\Delta v_r^{bp}$ conditioned on pairs with 
$r\in [r_c-\eta /2,r_c+\eta /2]$ at small 
$St$ closely follows the resulting expression for the relative velocity 
$\Delta v_r/u_\eta = St\,(\beta _b - \beta _p)a_{fr}/a_\eta$, where 
$a_{fr}$ is the average radial fluid acceleration at bubble and particle positions. We therefore conceptualise the resulting excess bp collision velocity as the local ‘turnstile mechanism’, since bubbles and particles respectively move in opposite directions, as illustrated in figure 10(c). An estimate for this effect can be obtained using 
$a_f \sim a_\eta$ in (2.9) and assuming a random orientation between the particle separation vector and the fluid acceleration, resulting in an average angle 
$\phi _{rdm} = 1\,\mathrm {rad} \approx 57^\circ$. This is shown by the dashed line in figure 10(a). While not strictly a prediction for 
$\Delta S_-$, this estimate is largely consistent with the data for 
$St<1$, suggesting that the enhanced approach velocity in this range is indeed due to the local turnstile effect. This is further consistent with the fact that for 
$St\geq 1$, the correlation of 
$\Delta v_r^{bp}$ with the local fluid acceleration is seen to vanish almost completely in figure 10(b). A positive 
$\Delta S_-$ is also predicted by the model of Zaichik et al. (Reference Zaichik, Simonin and Alipchenkov2010), but underestimates the value in the low-
$St$ regime. Beyond 
$St = 1$, 
$\Delta S_-$ flattens off at high values. Along with the observations made in figure 8(c), this suggests that the ‘path history’ effect active at these 
$St$ values is enhanced for the bp case. We can rationalise this in terms of a non-local ‘turnstile mechanism’ (see figure 10c), where differences in the particle velocities arise as bubbles and particles interact differently with larger and more energetic eddies in the flow. As a consequence, their local velocities differ more than for monodisperse cases at equal 
$St$. Both the local and non-local inertial effects discussed here enhance the relative approach velocity and therefore explain why 
$\varGamma ^{(STc)}_{bp} < \varGamma _{bp}$ as observed before.
(a) The excess bubble–particle effective radial collision velocity 
$\Delta S_- = S_-^{bp} - (S_-^{bb}+S_-^{pp})/2$ from simulations over a range of 
$St$. The factor 
$1/2$ for the dashed line accounts for the fact that 
$S_-$ is obtained by averaging only over the negative portion of 
$\mathrm {p.d.f.}(\Delta v_r|r_c)$ as defined in (1.3). (b) The radial component of the relative velocity conditioned on pairs with 
$r\in [r_c-\eta /2,r_c+\eta /2]$ binned by the local radial fluid acceleration at 
$Re_\lambda =175$. (c) Schematic illustrations of the local and non-local turnstile mechanisms.
4.2. Lagrangian statistics
To better examine the collision process, we adopt the Lagrangian point of view, where colliding pairs are tracked individually for separations 
$r \leq 100\eta$ up to collision. The corresponding statistics averaging over all pairs with the same 
$r$ are denoted by the suffix 
$\mathcal {L}$.
In figure 11, we show how 
$R^2|_\mathcal {L}$ and 
$\mathcal {S}^2|_\mathcal {L}$ vary as functions of 
$r$ for the colliding pairs. Consistent with the results in figure 5, variations are significantly stronger for 
$R^2|_\mathcal {L}$ (figure 11a) compared to 
$\mathcal {S}^2|_\mathcal {L}$ (figure 11b). In contrast to the monodisperse case, bubbles and particles need to leave their preferred flow regions in order to collide. While at 
$r/\eta \sim 100$ the unconditioned preferential concentration (indicated by the symbols in figure 11) is mostly recovered, a pronounced drift generally sets in for 
$r/\eta \lessapprox 50$. This might suggest that interaction with eddies at this intermediate range plays a role in bringing bubbles and particles together from their segregated locations. At 
$St = 0.1$, bubbles and particles appear to occupy similar flow regions for 
$r/\eta \lessapprox 10$ before colliding, i.e. in particular, the 
$R^2|_\mathcal {L}$ curve flattens at these scales, whereas the same is not the case at 
$St = 1$ and 
$St = 3$, where the collision distance is larger. These observations are consistent with the considerations regarding the local and non-local turnstile mechanisms made above. It is noteworthy that 
$R^2|_\mathcal {L}$ at 
$St \leq 1$ varies mostly for bubbles as they approach particles, while at 
$St = 3$, 
$R^2|_\mathcal {L}$ is almost constant for the bubbles but varies more for the particles. The latter case also stands out in the 
$\mathcal {S}^2|_\mathcal {L}$ plot as the slight decrease in strain values towards collision is not observed there.
Norm of (a) the rotation rates and (b) the strain rates of the flow at bubble and particle positions against pair separation conditioned on colliding pairs at 
$Re_\lambda = 175$ for 
$r \geq r_{c}$. Only the 
$St = 1$ case is shown for bubble–bubble and particle–particle pairs. Also plotted at 
$r\rightarrow +\infty$ is the unconditioned 
$R^2$ and 
$\mathcal {S}^2$ at the respective 
$St$.
In figure 12, we present velocity statistics conditioned on the colliding pairs. In particular, we consider the approach velocity 
$S_-|_\mathcal {L}$ (figure 12a) as well as the r.m.s. relative speed 
$(\Delta v)'|_\mathcal {L}$ (figure 12b). The approach velocities differ slightly in magnitude from the unconditioned results in figure 8, but the trends remain consistent. In particular, 
$S^{bp}_-|_\mathcal {L}$ is again largest for separations close to the respective collision distances for all 
$St$ shown. This is remarkable in view of the fact that 
$(\Delta v)'|_\mathcal {L}$ is actually largest for the bb case. This implies that 
$\Delta \boldsymbol {v}$ and 
$\boldsymbol {r}$ tend to be more anti-aligned for the bp case to yield the higher relative approach velocity.
(a) The effective approach velocity, and (b) the r.m.s. relative speed, conditioned on colliding pairs at 
$Re_\lambda = 175$.
To measure this, the angle between the separation vector and the relative velocity vector 
$\Delta \phi$ (see the inset of figure 13a) is used. For a head-on approach, 
$\Delta \phi =0^\circ$, while 
$\Delta \phi = 90^\circ$ means that the other particle is circling around in the rest frame of the collision partner. Note that as colliding particles must be approaching each other, 
$0^\circ \leq \Delta \phi <90^\circ$. As demonstrated in figure 13(a), the bubble–particle pairs are indeed the most anti-aligned, i.e. 
$\Delta \phi |_\mathcal {L}$ is lowest, relative to the other types of collisions near 
$r_{c}$. Figure 13(b) shows the ratio of the cosines of 
$\Delta \phi |_\mathcal {L}(r_c)$ as a measure for how much differences in alignment lead to a relative enhancement of 
$S_-(r_c)$. The ratio is up to 10 for bp relative to bb collisions, and up to 
${\sim }1.5$ compared to particle–particle collisions, and values for both cases are slightly higher at the higher 
$Re_\lambda$. When interpreting the smaller difference between pp and bp collisions here, it should be kept in mind that based on the discussion around figure 7, the collision mechanism differs between these two cases: whereas pp collisions appear predominantly shear-driven, where a good alignment may be expected, inertial effects, for which this is not necessarily the case, play a much more significant role for bp collisions.
(a) The angle between the separation and relative velocity vectors 
$\Delta \phi$ for 
$r \geq r_{c}$ at 
$Re_\lambda = 175$. Bubble–particle/bubble–bubble/particle–particle pairs are represented by green/blue/red lines as in figure 12. The inset shows the definition of 
$\Delta \phi$ in the rest frame of the particle on the right. (b) The extra angular contribution to the radial component of the bubble–particle collision velocity relative to the bubble–bubble and particle–particle cases.
4.3. Effect of lift, finite particle density and nonlinear drag
Our simulations are performed without the lift force using a finite particle density 
$\rho _p/\rho _f=5$ and a nonlinear drag law with 
$f_i = 1+0.169\,Re_i^{2/3}$. To test how sensitive our results are to these parameters, we ran addition simulations for 
$Re_\lambda =175$ changing one parameter at a time: adding the lift force 
$\boldsymbol {F_L} = -2/3\rho _f{\rm \pi} r_i^3(\boldsymbol {v}_i-\boldsymbol {u})\times (\boldsymbol {\nabla }\times \boldsymbol {u})$ to the right-hand side of (3.2), setting 
$\rho _p/\rho _f=\infty$, or using Stokes drag (
$f_i = 1$) for both phases. For all of these simulations, 
$r_c$ is chosen to be identical to the case when 
$\rho _p/\rho _f=5$. We compare the results for 
$\varGamma$ in figure 14, and those for 
$g(r_{c})$ and 
$S_-(r_{c})$ in figure 15. The plots show that the collision statistics are not sensitive to the lift force. Also more generally, the effect of the lift force is limited, with the most noticeable change occurring for the preferential sampling where we observed a slight decrease in 
$\tau _\eta \langle R^2\rangle _r$ at bubble clusters, consistent with Mazzitelli, Lohse & Toschi (Reference Mazzitelli, Lohse and Toschi2003). It is also evident that the particle density affects these quantities only marginally. Similarly, the results remain largely unchanged for the 
$f_i = 1$ case, with the exception of a slight increase in 
$\varGamma _{bp}$ due to a higher 
$S_-$ at the larger 
$St$. We furthermore note that other expressions for 
$f_i$ have been proposed in the literature (Schiller & Naumann Reference Schiller and Naumann1933; Clift, Grace & Weber Reference Clift, Grace and Weber2005; Wan et al. Reference Wan, Yi, Wang, Sun, Chen and Wang2020). These are compared in Appendix C, which shows that they do not differ significantly when 
$Re_i \lesssim 100$. Hence we conclude that our results would not be sensitive to whichever one of these parametrisations is chosen. Based on the findings in this section, it therefore appears appropriate to neglect lift, and to use the simplifications of linear drag and 
$\rho _p/\rho _f \to \infty$ for simulations and modelling approaches in the present parameter range.
The dimensionless collision kernel for different particle densities simulated, with the lift force included and for the case when the drag correction factor is 
$f_i = 1$.
(a) The RDF at collision distance, and (b) the effective radial collision velocity for different particle densities, with the lift force included and for the case when the drag correction factor is 
$f_i = 1$. The symbols follow figure 14.
5. Discussion and conclusion
We studied bubble–particle collisions in turbulence using a point-particle approach. Besides a critical appraisal of current models for the problem, our results highlight that the difference in the density ratios between bubbles and particles critically affects the collision statistics. An overview of the physical picture that emerges is presented in figure 16. For 
$St \to 0$, the shear mechanism is applicable, meaning that the collision velocity is well-described by the local fluid velocity gradient. Once 
$St$ reaches finite values, bubbles and particles concentrate preferentially in different regions of the flow, which leads to their segregation. This effect reduces the bubble–particle collision kernel 
$\varGamma _{bp}$ and is maximal for 
$St \sim 1$. The correction remains 
${O}(1)$, which is comparable to the enhancement of the collision kernel due to clustering for particle–particle cases, but much less than that for the bubble–bubble cases, respectively. For particle–particle collisions (evaluated at 
$r_c = r_b + r_p$), the relative approach velocities are found to be entirely consistent with the shear mechanism for the full range of 
$St$ explored here, such that the collision kernel in these cases is very well approximated by extending the shear-induced collision kernel to account for non-uniform particle concentration 
$\varGamma ^{(STc)}$. The same does not hold for bubble–particle and bubble–bubble collisions, where inertial effects additionally play a role and enhance the effective relative approach velocities beyond the shear scaling. This enhancement is strongest for the bubble–particle case, and we conceptualise this effect in terms of a ‘turnstile’ mechanism, which is related to the fact that fluid accelerations lead to opposing drift velocities relative to the fluid for bubbles and particles. This happens locally, i.e. on the same fluid element, if 
$St$ is sufficiently small for bubbles/particles to follow pathlines (
$St\leq 0.5$ based on the results in figure 10), but also appears to play a role at higher 
$St$, where the path history becomes more relevant. For 
$St \gtrapprox 2$, the ratio 
$\varGamma ^{(STc)}/\varGamma$ and the values of the effective radial approach velocity at contact 
$S_-(r_c)$ are comparable for bubble–particle and bubble–bubble collisions, which indicates that the ‘non-local’ turnstile is less prevalent there, in line with the fact that local particle/bubble velocities become increasingly random at larger 
$St$.
A sketch of the bubble–particle collision mechanisms with their approximate 
$St$ dependence.
Our method does not allow us to assess the high 
$St$ (kinetic-gas-like) regime since bubble sizes violate the point-particle assumption. Applying the Abrahamson (Reference Abrahamson1975) model pertinent to this regime to the present data vastly overpredicts the collision kernel, and the apparent improvement by adaptations of this framework is merely an artefact of using inconsistent velocity expressions. Our data for bubble–particle collisions are better approximated by the models of Saffman & Turner (Reference Saffman and Turner1956) and Zaichik et al. (Reference Zaichik, Simonin and Alipchenkov2010) when taking an additional correction for the inhomogeneous distribution into account. The predictions based on the theory of Zaichik et al. (Reference Zaichik, Simonin and Alipchenkov2010) generally overestimate 
$S_-$ and fail to reproduce the finding that 
$S_-^{bp}$ is largest in the present parameter range when comparing across species. These observations are not affected when using a linear drag relation, i.e. 
$f_i = 1$ in (3.2), or in the limit 
$\rho _p/\rho _f \to \infty$ as shown in § 4.3. We therefore could not reproduce the good agreement with the Zaichik et al. (Reference Zaichik, Simonin and Alipchenkov2010) theory for 
$S_-^{bp}$ reported in Fayed & Ragab (Reference Fayed and Ragab2013). Potentially, this is due to different choices for 
$St_b$, which we set to equal to 
$St_p$, whereas this parameter is kept constant in Fayed & Ragab (Reference Fayed and Ragab2013).
In our simulations, we employ the non-dimensional control parameters 
$St$ and 
$Re_\lambda$. In practice, the fluid and particle properties are controlled separately so using 
$r_b$, 
$r_p$ and 
$\varepsilon$ is more common. However, by definition, changing 
$\varepsilon$ affects both 
$St$ and 
$Re_\lambda$. To assess the overall effect, consider the region between 
$St = 0.1$ and 
$St = 0.5$, where the non-dimensional collision kernel 
$\varGamma _{bp}\tau _\eta /r_c^3$ decreases most sharply and drops by a factor of approximately 1.3. As the particle parameters are constant, 
$\tau _i = \textrm {const.}$, this is associated with a 5-fold decrease in 
$\tau _\eta$. According to the shear scaling, which appears to be applicable to our data, this increase in turbulence intensity implies also a 5-fold increase in 
$\varGamma _{bp}$. Hence there remains a net benefit from increasing the turbulence intensity even in this range, although the effect is reduced somewhat due to segregation.
Our results show that considering finite particle density, nonlinear drag and the lift force have limited effect on the bubble–particle collision rate. However, other relevant factors are yet to be explored. Most importantly, this concerns the buoyancy force and finite size effects in particular for the bubble motion and non-identical 
$St$ for bubbles and particles, as well as modelling the collision behaviour of bubbles and particles. More complex simulations as well as experiments are certainly needed to make precise predictions of the bubble–particle collision rate in realistic settings. However, as far as modelling approaches and a physical understanding are concerned, we believe that the most promising approach to disentangle the effect of additional factors (such as buoyancy) is to treat them as modifications to a simpler base case like the one studied here.
Acknowledgements
We thank R. Ostilla-Mónico and L. Jiang for useful discussions, as well as R. Yang, S. Yerragolam and C. Howland for help with visualisation and code compilation.
Funding
This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 950111, BU-PACT). The simulations are conducted with the Dutch National Supercomputers Cartesius and Snellius, and MareNostrum4 of the Barcelona Supercomputing Center.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
All data supporting this study are openly available from the 4TU.ResearchData repository at https://doi.org/10.4121/22147319.
Appendix A. Extension of the Kruis & Kusters (Reference Kruis and Kusters1997) model to particles with different densities
The original model by Kruis & Kusters (Reference Kruis and Kusters1997) is applicable to collisions of particles with arbitrary but equal density only. Here, we extend their framework to collisions of particles with different densities. In their model, Kruis & Kusters (Reference Kruis and Kusters1997) considered two limits: small and large 
$St$. We follow their approach to derive the corresponding expressions for particles with different densities. For convenience, we first rewrite (3.2) as
\begin{equation} \frac{{\rm d}\boldsymbol{v}_i}{{\rm d}t} = \gamma_i\,\frac{\mathrm{D}\boldsymbol{u}}{\mathrm{D}t} + \frac{\boldsymbol{u}-\boldsymbol{v}_i}{\tau_i}, \end{equation}
where 
$\gamma _i = {3\rho _f}/({2\rho _i + \rho _f})$ and 
$i = 1,2$ as in § 1.
A.1. Small $St$ limit
For small 
$St$, Kruis & Kusters (Reference Kruis and Kusters1997) based their model on Yuu (Reference Yuu1984), which considers a ‘local’ inertial effect due to the different response to fluid fluctuations of particles having different 
$St$ (originally termed the ‘accelerative mechanism’ in Kruis & Kusters (Reference Kruis and Kusters1997)) and the shear mechanism. The main contribution by Kruis & Kusters (Reference Kruis and Kusters1997) is employing a more accurate expression of the fluid Eulerian energy spectrum that also captures the dissipation range behaviour
\begin{equation} E_f(\omega) = \frac{4u'^2\xi}{2{\rm \pi}(\xi-1)}\left(\frac{T_{fL}}{1+T_{fL}^2\omega^2} - \frac{T_{fL}}{\xi^2 + T_{fL}^2\omega^2} \right), \end{equation}
where 
$\omega$ is the angular frequency,
\begin{equation} \xi = 2\left(\frac{T_{fL}}{\tau_L}\right)^2, \end{equation}
$T_{fL}=0.4L_f/u'$ is the fluid Lagrangian integral time scale, 
$L_f$ is the Eulerian longitudinal integral length scale, and 
$\tau _L^2=2u'^2/\langle ({\rm D}u/{\rm D}t)^2\rangle = 2u'^2/(1.16\varepsilon ^{3/2}\nu ^{-1/2})$ is the square of the Lagrangian time scale.
For the local inertial effect, the ensemble average of the square of the collision velocity is
\begin{equation} \langle[(\Delta v)^{(loc)}_{12}]^2\rangle = v_1'^2 + v_2'^2 - 2\langle v_1 v_2\rangle, \end{equation}
where 
$\langle \cdot \rangle$ denotes ensemble averaging, and
\begin{equation} v_i^2 = u'^2\,\frac{\xi}{\xi-1}\left[\frac{1 + \gamma_i^2 \tau_i/T_{fL}}{1 + \tau_i/T_{fL}} - \frac{1 + \gamma_i^2\xi\tau_i/T_{fL}}{\xi(1+\xi\tau_i/T_{fL})}\right] \end{equation}as derived in Kruis & Kusters (Reference Kruis and Kusters1997). The integral giving the particle velocity correlation term follows Yuu (Reference Yuu1984), whose model has been extended to particles with unequal densities in Ngo-Cong et al. (Reference Ngo-Cong, Nguyen and Tran-Cong2018). The integral in Ngo-Cong et al. (Reference Ngo-Cong, Nguyen and Tran-Cong2018) reads
\begin{align} \langle v_1 v_2\rangle = \int_{0}^{+\infty} \frac{(\tau_1^{{-}2}+\gamma_1 \omega ^2) (\tau_2^{{-}2}+\gamma_2 \omega ^2)+\tau_1^{{-}1} \tau_2^{{-}1} (\gamma_1-1) (\gamma_2-1) \omega^2}{(\tau_1^{{-}2}+\omega ^2) (\tau_2^{{-}2}+\omega^2)}\,E_f(\omega) \,{\rm d}\omega. \end{align}Note that although Kruis & Kusters (Reference Kruis and Kusters1997) found a misprint in the integral given by Yuu (Reference Yuu1984), this error has not propagated to Ngo-Cong et al. (Reference Ngo-Cong, Nguyen and Tran-Cong2018). Performing the integration yields
$$\begin{gather} \frac{\langle v_1 v_2\rangle}{u'^2} = 1 + \frac{\xi}{(\xi-1)(\tau_1^{{-}1} + \tau_2^{{-}1})}\left[\frac{(\gamma_1-1)(\tau_2^{{-}1}+\tau_1^{{-}1}\gamma_2)}{1+\tau_1^{{-}1}T_{fL}} + \frac{(\tau_1^{{-}1} + \tau_2^{{-}1}\gamma_1)(\gamma_2-1)}{1+\tau_2^{{-}1}T_{fL}}\right. \nonumber\\ \left.{}- \frac{(\gamma_1-1)(\tau_2^{{-}1} + \tau_1^{{-}1}\gamma_2)}{\tau_1^{{-}1}T_{fL} + \xi} - \frac{(\tau_1^{{-}1} + \tau_2^{{-}1}\gamma_1)(\gamma_2 - 1)}{\tau_2^{{-}1}T_{fL} + \xi} \right], \end{gather}$$
which is different from the corresponding expression in Ngo-Cong et al. (Reference Ngo-Cong, Nguyen and Tran-Cong2018) solely because of the usage of another form of 
$E_f(\omega )$ as given by (A2) in (A6).
For the collision velocity due to the shear mechanism, Kruis & Kusters (Reference Kruis and Kusters1997) used the result in Yuu (Reference Yuu1984):
\begin{equation} \langle[(\Delta v)^{(shr)}_{12}]^2\rangle = \frac{\varepsilon}{5\nu}\left(\frac{v_1'^2}{u'^2}\,r_1^2 + \frac{v_2'^2}{u'^2}\,r_2^2 + 2\,\frac{\langle v_1 v_2 \rangle}{u'^2}\,r_1r_2\right). \end{equation}
Using (1.5) and allowing 
$\gamma _1 \neq \gamma _2$ gives
\begin{equation} \langle[(\Delta v)^{(shr)}_{12}]^2\rangle = \frac{3\varepsilon}{5}\left(\frac{v_1'^2}{u'^2}\,\frac{\tau_1}{\gamma_1} + \frac{v_2'^2}{u'^2}\,\frac{\tau_2}{\gamma_2} + 2\,\frac{\langle v_1 v_2 \rangle}{u'^2}\sqrt{\frac{\tau_1\tau_2}{\gamma_1\gamma_2}}\right). \end{equation}The overall collision kernel is then given by
\begin{equation} \varGamma_{12} = \sqrt{8{\rm \pi}}\,r_c^2\,\sqrt{\langle[(\Delta v)^{(loc)}_{12}]^2\rangle + \langle[(\Delta v)^{(shr)}_{12}]^2\rangle}. \end{equation}
Note that in contrast to the expression in Kruis & Kusters (Reference Kruis and Kusters1997), (A10) does not carry a factor 
$1/\sqrt {3}$ as 
$u'$ is defined as the single-component r.m.s. fluid velocity.
A.2. Large $St$ limit
The Kruis & Kusters (Reference Kruis and Kusters1997) expression of 
$\varGamma$ for large 
$St$ is essentially based on Williams & Crane (Reference Williams and Crane1983) except that the added mass term is retained. As in § A.1, the densities of the different species 
$\gamma _i$ are assumed to be distinct. In this limit, the shear mechanism is not considered, so
\begin{equation} \langle(\Delta v)^2_{12}\rangle = v_1'^2 + v_2'^2 - 2\langle v_1 v_2\rangle \end{equation}and
\begin{equation} \varGamma_{12} = \sqrt{8{\rm \pi}}\,r_c^2\,\sqrt{\langle(\Delta v)^2_{12}\rangle}. \end{equation}
Furthermore, a simpler version of 
$E_f(\omega )$, namely
\begin{equation} E_f(\omega) = \frac{4u'^2}{2{\rm \pi}}\left(\frac{T_{fL}}{1+T_{fL}^2\omega^2}\right) \end{equation}
is used, which corresponds to an exponentially decaying fluid velocity autocorrelation function ‘moving with the mean flow’ 
$R_E^*(t) = \exp ({-|t|/T_{fL}})$ that does not account for the dissipation range. Kruis & Kusters (Reference Kruis and Kusters1997) showed that
\begin{equation} v_i'^2 = u'^2\,\frac{1 + \gamma_i^2\tau_i/T_{fL}}{1 + \tau_i/T_{fL}}, \end{equation}which leaves only the particle velocity correlation to be determined.
To determine 
$\langle v_1v_2\rangle$, we integrate (A1) for 
$t\in (-\infty,0]$ and exercise the freedom that 
$t = 0$ can be chosen at any instant, to obtain
\begin{equation} v_{ij}(\boldsymbol{x_i},t) = \gamma_i\,u_{ij}(\boldsymbol{x_i},t) + \frac{1-\gamma_i}{\tau_i}\int_{0}^{+\infty}u_{ij}(\boldsymbol{x_i},t-\phi) \exp\left({-\frac{\phi}{\tau_i}}\right){\rm d}\phi, \end{equation}
where 
$j = x,y,z$ are the individual components of the corresponding vector, and 
$\boldsymbol {x_i}$ is the position vector. Then using the approximation (Williams & Crane Reference Williams and Crane1983)
\begin{equation} \langle u(\boldsymbol{x_1},t')\,u(\boldsymbol{x_2},t'') \rangle \approx R_E^*(t'-t'')\int_{0}^{+\infty}E_f(\omega)\cos\left[\omega\,\widehat{\Delta v}\left(t - \frac{t'+t''}{2}\right)\right]{\rm d}\omega, \end{equation}
where 
$\widehat {\Delta v}$ is a non-dimensional particle relative velocity, we have
\begin{align} &\langle v_1 v_2 \rangle =
\gamma_1\gamma_2u'^2 +
\frac{(1-\gamma_1)(1-\gamma_2)}{\tau_1\tau_2}
\int_{0}^{\infty}\int_{0}^{\infty}\int_{0}^{\infty}
\exp\left({-\frac{|\psi-\phi|}{T_{fL}}}\right)\frac{2}{\rm \pi}\,u'^2\,\frac{T_{fL}}{1+\omega^2T_{fL}^2}
\nonumber\\ &\quad \times\cos{\left[\omega\,\widehat{\Delta
v}\left(\frac{\psi+\phi}{2}\right) \right]}
\exp\left({-\frac{\psi}{\tau_1} -
\frac{\phi}{\tau_2}}\right){\rm d}\omega \,{\rm
d}\psi \,{\rm d}\phi \nonumber\\ &\quad +
\frac{(1-\gamma_1)\gamma_2}{\tau_1}\int_{0}^{\infty}\int_{0}^{\infty}
\exp\left({-\frac{|\psi|}{T_{fL}}}\right)\frac{2}{\rm \pi}\,u'^2\,\frac{T_{fL}}{1+\omega^2T_{fL}^2}\cos\left({\omega}\,\widehat{\Delta
v}\,\frac{\psi}{2}\right)
\exp\left({-\frac{\psi}{\tau_1}}\right){\rm
d}\omega \,{\rm d}\psi \nonumber\\ &\quad +
\frac{\gamma_1(1-\gamma_2)}{\tau_2}\int_{0}^{\infty}\int_{0}^{\infty}
\exp\left({-\frac{|\phi|}{T_{fL}}}\right)\frac{2}{\rm \pi}\,u'^2\,\frac{T_{fL}}{1+\omega^2T_{fL}^2}\cos\left({\omega}\,\widehat{\Delta
v}\,\frac{\phi}{2}\right)
\exp\left({-\frac{\phi}{\tau_2}}\right){\rm
d}\omega \,{\rm d}\phi.
\end{align}
As 
$\tau _1,\tau _2 \gg 1$, 
$\exp (-\psi /\tau _1 - \phi /\tau _2) \approx 0$ except when 
$\psi = \phi$. Hence we retain this case only in the arguments of cosine and 
$\exp (-\psi /\tau _1 - \phi /\tau _2)$. Furthermore, we replace 
$\widehat {\Delta v}$ by 
$\widehat {\Delta v_X} = \sqrt {(v_1'^2 + v_2'^2)/u'^2}$ as the particle velocity correlation should be weak at large 
$St$. Thus
\begin{align} &\frac{\langle v_1 v_2
\rangle}{u'^2} = \gamma_1\gamma_2 +
\frac{(1-\gamma_1)(1-\gamma_2)}{\tau_1\tau_2}
\int_{0}^{\infty}\int_{0}^{\infty}\int_{0}^{\infty}
\exp\left({-\frac{|\psi-\phi|}{T_{fL}}}\right)
{\rm
d}\psi\,\frac{2T_{fL}}{{\rm \pi}(1+\omega^2T_{fL}^2)}
\nonumber\\ &\quad \times\cos{(\omega\,\widehat{\Delta v_X}\,\phi)}
\exp\left[{-\phi\left(\frac{1}{\tau_1} +
\frac{1}{\tau_2}\right)}\right]{\rm
d}\omega \,{\rm d}\phi \nonumber\\ &\quad +
\frac{(1-\gamma_1)\gamma_2}{\tau_1}\int_{0}^{\infty}\int_{0}^{\infty}
\exp\left({-\frac{|\psi|}{T_{fL}}}\right)\frac{2T_{fL}}{{\rm \pi}(1+\omega^2T_{fL}^2)}\cos\left(\omega\,\widehat{\Delta
v_X}\,\frac{\psi}{2}\right)
\exp\left({-\frac{\psi}{\tau_1}}\right) {\rm
d}\omega \,{\rm d}\psi \nonumber\\ &\quad +
\frac{\gamma_1(1-\gamma_2)}{\tau_2}\int_{0}^{\infty}\int_{0}^{\infty}
\exp\left({-\frac{|\phi|}{T_{fL}}}\right)\frac{2T_{fL}}{{\rm \pi}(1+\omega^2T_{fL}^2)}\cos\left(\omega\,
\widehat{\Delta v_X}\,\frac{\phi}{2}\right)
\exp\left({-\frac{\phi}{\tau_2}}\right) {\rm
d}\omega \,{\rm d}\phi.
\end{align}
Integrating first over 
$\omega$ and then the remaining variables finally results in
\begin{align} \frac{\langle v_1 v_2 \rangle}{u'^2} & = T_{fL}\,\frac{(1-\gamma_1)(1-\gamma_2)}{\tau_1\tau_2} \left[\frac{2}{\tau_1^{{-}1} + \tau_2^{{-}1} + \widehat{\Delta v_X}\,T_{fL}^{{-}1}} - \left(\frac{1}{\tau_1} + \frac{1}{\tau_2} + \frac{\widehat{\Delta v_X} + 1}{T_{fL}} \right)^{{-}1} \right] \nonumber\\ &\quad + \frac{(1-\gamma_1)\gamma_2}{\tau_1}\left(\frac{1}{T_{fL}} + \frac{1}{\tau_1} + \frac{\widehat{\Delta v_X}}{2T_{fL}} \right)^{{-}1} + \frac{\gamma_1(1-\gamma_2)}{\tau_2}\left(\frac{1}{T_{fL}} + \frac{1}{\tau_2} + \frac{\widehat{\Delta v_X}}{2T_{fL}} \right)^{{-}1} \nonumber\\ &\quad + \gamma_1\gamma_2. \end{align}Appendix B. Verification of point-particle code
To verify our code for the suspended phase, we compared the results of the infinitely heavy particle cases to the literature. As demonstrated by figure 17, the RDF 
$g(r)$, the variance of the radial component of the relative velocity 
$S_{2\parallel }^{pp}(r)$ as defined by (4.1), and the collision kernel 
$\varGamma$ agree very well with data from Ireland et al. (Reference Ireland, Bragg and Collins2016) and Voßkuhle et al. (Reference Voβkuhle, Pumir, Lévêque and Wilkinson2014), especially considering the difference in 
$Re_\lambda$ compared to the literature data in figures 17(a,b).
Particle–particle collision statistics of infinitely heavy particles: (a) the RDF 
$g_{pp}(r)$ at 
$St = 1$; (b) the variance of the radial component of the relative velocity 
$S_{2\parallel }^{pp}(r)$ at 
$St = 1$; and (c) the collision kernel 
$\varGamma _{pp}$ when taking 
$\rho _p/\rho _f=250$ and 
$r_c = 2r_p$.
Appendix C. Comparison of different drag parametrisations
Various expressions of the drag correction factor 
$f_i$ have been proposed in the literature. Figure 18 shows that they are mostly similar for 
$Re_i \lessapprox 100$.
Various drag correction factors 
$f_i$ from the literature compared to the one by Nguyen & Schulze (Reference Nguyen and Schulze2004) used in this study (solid green line).
Open access
