## 2 Context and Summary of Numerical Simulations  
The simplest and most extensively studied type of turbulence is statistically homogeneous isotropic turbulence of an incompressible fluid. In such configuration, it is well known that small-scale fields of velocity can be thought of as being more or less homogeneous and isotropic. Small-scale vortex structures are responsible for dissipation of turbulent energy and play a key role in the accumulation (clustering) of particles in a turbulent flow. Study of isotropic turbulence is thus of fundamental importance for both single-phase and two-phase flows and, therefore, it has been the configuration adopted in the present work. In this work, the simulation of the turbulent behaviour of the flow has been performed by the Lattice Boltzmann Method with the in-house code described in Reference [38] and also employed in Reference [38] to study the effect of inter-particle collisions on several clustering measures of inertial particles in turbulent flows. For the details of the implementation of the Lattice Boltzmann Method (LBM) forcing scheme for generating the homogeneous isotropic turbulence and collision detection algorithm, the reader is referred to Reference [38]. Therefore, only an abridged summary of the main flow characteristics is mentioned here. We focus on the dynamics of colliding particles in an isotropic homogeneous flow in a triple-periodic box of size \(L=0.128\) m. The fluid domain is uniformly discretised by \(128^{3}\) grid points. Forcing parameters have been adjusted to guarantee a good resolution of all flow scales, providing a grid space slightly smaller than the Kolmogorov length scale \(\eta_{K}\). Additionally, particle size is much smaller than \(\eta_{K}\) so they are treated as point particles; moreover, the simulated particle volume fractions are sufficiently small to keep the two-phase flow dilute, so the back-reaction force of the particles on the flow is neglected. For the same reason, only binary elastic collisions are considered in the framework of the hard sphere approach. A brief summary of the employed collision model and detection algorithm is presented in Appendix A. Particle trajectories are computed by the simplified particle equation of motion, retaining only the contribution of drag term, in combination with the Schiller and Naumann correlation [39; 40]. Domain discretisation and particle volume fraction values have been chosen to keep the computational effort limited as the direct detection of particle collisions is very costly. The main parameters employed in the fluid simulation are presented in Table 1.  
As an illustration, Figure 1 shows the non-dimensional three-dimensional (3D) energy spectrum \(\varepsilon_{x}\) obtained in the present DNS. For comparison, the experimental results of [PERSON] and Corrsin [41] and the Kolmogorov spectrum are also presented. The considered particles have a fixed diameter, \(d_{p}=200\)\(\upmu\)m, which is more than five times smaller than the Kolmogorov length scale \(\eta_{K}\); therefore, particle size is well in the dissipative range. Three different average particle volume fractions, \(\alpha_{p}=0.01\%\), \(0.05\%\) and \(0.1\%\), have been simulated. Particle volume fraction is varied by changing the number of tracked particles: for \(\alpha_{p}=0.01\%\) there are 50,066 particles, for \(\alpha_{p}=0.05\%\) there are 250,329 particles and for \(\alpha_{p}=0.1\%\), there are 500,658 particles. Particle properties are summarised in Table 2.  
Inertial particles' dynamics depend on Stokes number, which characterises particle inertia and is equal to the ratio between the average response time of particles, \(\left<\tau_{p}\right>\), and the Kolmogorov time microscale, \(\tau_{K}\), \(St_{K}=\left<\tau_{p}\right>/\tau_{K}\).