Lagrangian statistics for fluid particles and bubbles in turbulence

The dispersion of bubbles in homogeneous and isotropic turbulence is numerically examined. The fluid velocity field evolves according to the Navier–Stokes equations that are solved by direct numerical simulation. The bubble paths are followed by Lagrangian tracking. The aim of the work is to quantify dispersion properties of bubbles in a regime ranging from low- to high-turbulence velocity fluctuations as compared to the bubble velocity scale, and to compare them to those of fluid particles. Moreover, the forces that are relevant for the bubble dispersion are analysed and their probability density functions, as well as their intermittency characteristics, are compared to those of fluid particle accelerations.


Introduction
An important aspect of turbulent flows is their ability to enhance the transport and the mixing of scalar quantities suspended within them. The dispersion of particles and bubbles in homogeneous and isotropic turbulence is of particular interest.
The picture that can be drawn from both experimental and numerical studies [9]- [12] is that fluid particles, in general, disperse more efficiently than heavy particles or bubbles whose diffusion is reduced by the drift along the gravity direction. Yet, a full understanding of the forces that are relevant in the phenomenon is lacking. For instance, whereas strong gravity effects are known to reduce heavy particle diffusion as compared to fluid particles, the role played by inertia is less evident, even though it is often found to increase the diffusivities [11]- [13].
In the present work, we focus on air microbubbles of diameter d ∼ 10 −4 m rising in uncontamined water (the fluid viscosity is ν = 10 −6 m 2 s −1 and the gravity is g = 9.8 m s −2 ). We numerically examine the bubble motion and dispersion in homogeneous and isotropic turbulence. Particular emphasis is focused on the similarities and differences between bubbles and fluid particles. Direct numerical simulation is employed to simulate the fluid velocity field and the Lagrangian tracking of the trajectories is used as an efficient tool for probing the statistics.
The paper is organized as follows: in section 2 the numerical method, extensively described in a former paper [8], is briefly reported. In section 3 the theory of diffusion of fluid particles and bubbles is reviewed and compared to the results of our study. Relative diffusion is addressed in section 4. The statistical characteristics of the fluid acceleration along particle trajectories and of the forces acting on bubbles are assessed in section 5. The high degree of intermittency revealed by their probability density functions is further investigated in section 6. Here, Lagrangian velocity structure functions, both for fluid particles and for bubbles, and in particular their scaling properties are addressed. Section 7 contains conclusions.

Flow simulation
The computational domain is a cube of side L 0 = 2π, consisting of N 3 = 128 3 mesh points with periodic boundary conditions. The flow velocity u(x, t) evolves according to the incompressible Navier-Stokes equation: where p(x, t) represents the pressure, f L (x, t) a large scale forcing and f b (x, t) the small scale bubble forcing. Equation (1) is solved in wavenumber k space, by using the pseudo-spectral method [8]. All flow scales from the integral, where the energy is introduced in the system, down to the Kolmogorov scale η, where it is dissipated by viscous effects, are resolved. A statistically stationary state is substained by the forcing (in wavenumber space), f L (k, t) = u(k, t) k∈K in |u i (k, t)| 2 , k ∈ K in (2) and f L (k, t) = 0 otherwise. f L (k, t) and u(k, t) are, respectively, the Fourier transform of f L (x, t) and u(x, t). K in is a subset of small wavenumbers defined according to: K in = {k such that(L 0 /(2π))k = ±(−1, 2, 2), ±(2, −1, −1) + permutations}. This type of forcing has been successfully applied in many of our previous simulations [9], [14]- [16]. Units are fixed by k min = 1 and = 1 for f b (x, t) = 0. Then good resolution down to the Kolmogorov scale is achieved when ν = 0.007. The Kolmogorov length, time, and velocity scales can be evaluated according to: η = (ν 3 / ) 1/4 , τ k = (ν/ ) 1/2 and v k = ( ν) 1/4 , respectively. Other relevant quantities are the large scale r.m.s. fluid velocity fluctuations u 0 , the large Eddy turnover time τ 0 , the integral length scale L 11 , 1 the Taylor scale λ and the Taylor Reynolds number Re λ .
(∂ x u x ) 2 1/2 and Re λ = u 0 λ ν , 1 The integral scale L 11 reported in table 1 (which is the length at which the Eulerian velocity autocorrelation function decreases to 1/e of its initial value) is larger with respect to the scale computed according to L 11 = 1 u 2 x π r x =0 u x (0)u x (r x ) r x , because the velocity autocorrelation function u x (0)u x (r x ) in turbulence simulations driven by the forcing (2) presents a negative loop on the large scales: for r x = π it is u x (0)u x (r x ) / u x (0) 2 ≈ −0.17, whereas for larger r x the correlation increases again due to the periodic boundary conditions.
The statistically stationary value of Re λ in the present simulation is 62. A list of the fluid parameters can be found in table 1. Two different cases are investigated in this paper: (i) the one-way coupling case, in which bubbles do not force the flow, i.e., f b (x, t) = 0 in equation (1), and (ii) the two-way coupling case, in which each bubble behaves as a point force [8] on the carrier flow. Where N b is the total number of bubbles, y i (t) is the instantaneous position of bubble i, g is the gravity, directed along the negative z-axis and δ is the Dirac δ-function [8], [17]- [19]. All bubbles are identical spheres of radius a and volume V b = 4πa 3 /3.

Bubble motion equation
The equation of motion for a small spherical bubble, immersed in the fluid field u, is modelled as where v(t) stands for the bubble velocity, y(t) the bubble's instantaneous position, ω = ∇ × u the fluid vorticity and τ b the bubble response time that is connected with the bubble-rise velocity v T in still fluid by v T = 2gτ b , with g = |g|. The fluid material derivative is evaluated at the bubble position. When the bubble Reynolds number Re = 2|u − v|a/ν is less than 1, the bubble timescale is given by τ b = a 2 /(6ν) [20,21]. The physical meaning of the four terms on the right hand side of equation (4) is, from left to right, as follows: added mass, drag, buoyancy and lift force. For a detailed discussion of these effective forces we refer to [22]- [27]. In [27] we have also discussed the generalization of equation (4) to a spherical particle of any mass, and in particular to heavy particles. Note in particular that for the lift force, a lift coefficient of C L = 1/2 has been assumed. This is an approximation valid only in the large bubble Reynolds number limit [28] which has been well confirmed through numerical simulations [29]. For smaller bubble Reynolds number, the lift coefficient depends both on the local vorticity, the local shear and the bubble Reynolds number itself. In pure vortical flow and for bubble Reynolds numbers smaller than ∼5, the lift coefficient can even become negative [23,24]. In pure shear flow, on the other hand, the lift coefficient remains positive. In this paper, we ignore this complication and stick to a constant lift coefficient C L = 1/2. This implies that the exact numbers resulting from our numerics should be taken with some care. The focus of this work is the mechanisms and the trends, not the exact numeral values and we feel that taking C L = 1/2 is still a better approximation than ignoring the lift force altogether (C L = 0) as done in many previous studies.

Diffusion in homogeneous and isotropic turbulence
The dispersion tensor of fluid particles is defined by where x i (x 0 , t) is the ith component of the position vector at time t of a fluid particle that was at x 0 at the initial time t = t 0 , x i0 is the ith component of the initial position, and the brackets · indicate ensemble average (over all particle locations). We review here in short the analysis of the time evolution of the dispersion tensor [1,3]. In the following, v(x 0 , t) indicates the Lagrangian velocity of the particle. It equals the Eulerian velocity field at the fluid particle location, v(x 0 , t) = u(x(x 0 , t), t). The particle equation of motion is Integration gives Therefore, the diagonal elements of the dispersion tensor are Here and in the following, repeated indices are not considered to be summed. When the turbulence is stationary and homogeneous, the ensemble average can be moved inside the integral and substituted with the time average. Moreover, the limits of integration can be changed, by recognizing that the integration plane is a square. Thus and, by changing the variables (t , t ) The time average velocity correlation depends only on the time difference τ. We can express it in terms of the Lagrangian velocity autocorrelation coefficient, defined by Note that R i L (τ) is symmetric under the transformation τ → −τ. Therefore, From equation (12), the asymptotic results for both small and large time intervals t − t 0 can be easily obtained. When t → t 0 , then R i L (τ) 1 for all τ, and For large time intervals where T i L is the Lagrangian integral timescale in the ith direction. As a consequence that is, the diffusion in the three directions x, y and z increases linearly with time. It can be quantified by the diffusivity Due to isotropy the diffusivity is the same in all directions and the index i can be dropped in the last equation. Therefore, in the following, we will calculate the average value of the fluid diffusivity along x, y, and z and indicate it with T L is the Lagrangian integral timescale and the fluid r.m.s. velocity fluctuations u 0 enter in the equation because of the equality between Lagrangian time average and Eulerian space average. Similar reasoning as for fluid particles can be carried out for heavy particles or for gas bubbles, provided that one studies the dispersion with respect to the average position x(x 0 , t) [5], as heavy particles/bubbles sink/rise due to gravity. In the present case x(x 0 , t) ∝ (0, 0, v T (t − t 0 )), with v T the bubble-rise velocity in quiescent water. The prefactor is smaller than 1 as the turbulent fluctuations reduce the average bubble-rise velocity [8,9]. Thus, defining and indicating with X i (t) its ith component, the bubble dispersion tensor is The prediction for the behaviour of the diagonal components of the tensor at small and large time intervals (t − t 0 ) is, again, respectively quadratic and linear in (t − t 0 ) (see equations (13) and (15)). Dispersion properties can be quantified in three different ways, namely through (i) the diffusivity (ii) the ensemble average bubble velocity fluctuations Here R i L (τ) is the bubble Lagrangian velocity autocorrelation coefficient T b,i L is related to the diffusivity by In a bubbly flow there is no isotropy due to buoyancy. Therefore, we separately compute the lateral diffusivity D b xy = (D b x + D b y )/2 and the one in the direction of gravity, D b z . The values measured for bubbles will be compared to the ones for fluid particles.
Bubble dispersion has been approached analytically by Spelt and Biesheuvel [6,7] (henceforth indicated as SB97 [6] and SB98 [7]) and the predictions have been compared to the numerical results obtained by kinematic simulations of turbulence. The various regimes have been characterized by two dimensionless parameters, namely that represents the ratio of the r.m.s. fluid velocity fluctuations u 0 to the bubble-rise velocity v T , and with L 11 the integral scale of the turbulence and τ b the bubble response time. Large µ values indicate that the bubbles rapidly adapt their velocities to the one of the flow.
In the first paper (SB97), large bubble-rise velocities (i.e., β 1) have been considered. Moreover, u 0 /v T µ v T /u 0 was assumed. In the second paper (SB98), the case of large scale turbulence has been addressed (i.e., µ 1 and β 1). Analytical expressions for the diffusivities were estimated in the limit µ → ∞ and also when the first order correction in 1/µ was taken into account.
The range of parameters analysed in the present study corresponds to the regime investigated in the second work (SB98), so that the comparison to the analytical theory is straightforward. Furthermore, the advantage of simulating turbulence starting from the Navier-Stokes equations (as compared to kinematic simulations [7]) is to have a realistic description of both the small and the large flow scales. Therefore, it is possible to extend the analysis to a wide range of scales, in which the bubbles respond to both small and large scale turbulence fluctuations. The kinematic simulations employed in the previous papers do not properly reproduce the small scale structures of the vorticity field; as a consequence the approach is suitable only to small β values, i.e., to bubbles moving fast across the vortices, with only little interaction.
Small scale fluctuations are particularly important for the Lagrangian and Eulerian velocity autocorrelation at small time lags. We will now compare these two types of autocorrelations. The first one has already been defined in equation (11). The latter one is defined by where r is the spatial separation vector and the average is over space and time. We define . According to Tennekes [30], at small τ, owing to the advection or sweeping of the small eddies by the large scales, the Lagrangian autocorrelation function lies above the Eulerian one. However, it was shown by numerical simulation [12] that this effect may be masked at small Re, whereas it is certainly detectable at large Re. In our simulation, with Re λ = 62, we can barely see any crossover between the autocorrelation functions, as shown in figure 1. In kinematic simulations, the sweeping effect of the small by the large eddies is not modelled, and no crossover exists independent of the Re.
The Eulerian (τ E ) and Lagrangian (τ L ) time microscales can be computed by evaluating the curvature at time τ = 0 of the corresponding autocorrelation functions. This measure is very sensitive to numerical errors. We estimate where the velocities u i , v i and the time derivatives ∂u i /∂t, dv i /dt are evaluated, respectively, at fixed position in space for the Eulerian timescale and along the fluid particle trajectory for the Lagrangian one. In the present simulation, after averaging the results along the three directions, it is τ E = 0.33 and τ L = 0.35, thus τ L /τ E 1.06, which is in agreement with previous numerical results [32], at comparable Re λ .

Fluid particle diffusion
In order to fix reference values for the analysis of bubbles, we measure the dispersion of fluid particles by tracking the motion of N p = 40 960 particles that evolve according to equation (6). The velocity at the particle position is evaluated by third-order Taylor series-based interpolation [32] and the equation is advanced in time by explicit forward Euler scheme. The Courant number is C = t(3u 2 0 /2) 1/2 / x 1/30, with t the time step of the simulation and x = L 0 /N the mesh size. The code is run for several large-eddy turnover times until the values of large scale quantities, such as the total fluid energy and the viscous dissipation, become statistically stationary. Then all transient effects have decayed. After that particle statistics is accumulated for about five large-eddy turnover times.
The fluid diffusivity (see equation (17)) calculated in the present simulation is D f = 0.71u 0 L 11 . The Lagrangian integral timescale is thus: The Eulerian integral time can be estimated as T E = L 11 /u 0 ≈ 0.75 and the ratio T L /T E ≈ 0.7.

Bubble diffusion
We focus now on bubble diffusion. Several simulations are run with bubbles in the same carrier flow, i.e., the initial values of u 0 and L 11 are equal for all simulations. Later they can be modified by the two-way coupling effects. The trajectories and velocities of N b bubbles are tracked. The bubble response time is fixed (τ b = τ k /10, implying a/η ≈ 0.8) in all simulations. Each run is characterized by a different value of the bubble-rise velocity, chosen in the range v k v T 8v k . Consequently, each simulation corresponds to a different pair of parameters β = u 0 /v T and µ = L 11 /τ b v T . The smallest value of µ reached is 46, thus the first-order approximation (µ → ∞) of the analytical derivation presented in SB98 is likely to hold.
Both the one-way and two-way coupling regimes are investigated. In the first case, the bubbles move in the flow according to their equation of motion (4), without interacting, i.e., f b (x, t) = 0 in equation (1). The total number of bubbles tracked is N b = 40 960. In the twoway coupling regime, the bubble's feedback on the flow, as formulated in equation (3), is taken into account. In the latter simulations, the motion of N b = 144 000 bubbles is tracked, so that the void fraction is α = 1.6%. The regimes addressed correspond to microbubbles of diameters d ≈ 120-250 µm in water (ν = 10 −6 m 2 s −1 and g = 9.8 m s −2 ). All observables are measured separately along the three directions, x, y, and z. Then the values along x and y are averaged in order to obtain the lateral observable that is different from the vertical one (along z).
In figure 2, the vertical and lateral diffusivities of the bubbles are plotted as functions of β. Both of them are smaller than the diffusivity of fluid particles that is indicated by the broken straight line in the figure. We find the largest deviation for the smallest β value, where the bubblerise velocity is largest. The explanation can be sought in the so-called 'crossing trajectory effect' [33,34]: bubble trajectories 'cross' the trajectories of fluid particles, because of the gravity, and therefore rapidly leave those flow regions in which the velocity is highly correlated. The consequence is that the bubble velocity quickly loses the memory of its previous values by sampling flow velocities that are more and more decorrelated. Equations (21)- (23) show that the loss of velocity correlation reduces the diffusivity. The faster the bubbles rise (small β), the faster the decorrelation occurs. Furthermore, according to an effect related to the fluid incompressibility, the lateral diffusivity is smaller than the vertical [12].
The qualitative trends of the vertical and lateral diffusivities are in good agreement with the analytical results, see figure 1 of SB98. However, if we compare the results quantitatively, we find that the ratio of the lateral diffusivity to the fluid particle diffusivity is much smaller than theoretically predicted by SB98. This result was already obtained by the authors themselves, when studying bubble motion in turbulence kinematic simulations.
In order to explain the discrepancy between numerics and theory, we have a close look at the hypotheses made in order to derive formulae for the diffusivities. First, it is assumed that the probability density function (pdf) for the bubble displacement X(t) (see equation (18)) is Gaussian. In figure 3, we present the pdf(X(t)) evaluated at time t = 1.07L 11 /u 0 . The function is normalized in order to have unitary variance. For comparison, a Gaussian function with the same mean and variance is also plotted. The figure indicates that, although there are some deviations, the Gaussian approximation is reasonable. The second hypothesis of the analytical approach is Corrsin's independence conjecture [35], according to which particles experience unbiased fluid velocities, e.g. that are randomly distributed. This is clearly not the case in turbulent bubbly flows, where bubbles locally accumulate in eddies, owing to inertial forces, and in particular on their downflow side, mainly owing to the lift force [8,9,17]. The deviation from the analytical theory can thus be attributed to the relevance of inertial forces in bubble dynamics. Yet, the role played by vortex trapping and by the preferential sweeping in downward fluid velocity regions has to be still clarified.
To this end we focus on the diffusivity for β ∼ 0.7. In table 2 results are presented for three cases: (i) τ b = τ k /10, with lift force in the bubble equation of motion (4), (ii) τ b = τ k /10, without lift force and (iii) τ b = τ k , with lift force. For cases (i) to (iii) it is v T = 6v k . For completeness the table also contains the ratio of the enstrophy = 1 2 ω · ω averaged over all bubble locations normalized with the average flow enstrophy b / , the fraction N − /N b of the total number of bubbles that are located in downflow regions, and the relative reduction of the bubble-rise velocity in turbulence in comparison to quiescent fluid, ( v z b − v T )/v T . These observables indicate that, even if the turbulence intensity is rather low (i.e., β < 1), the main characteristics of bubble motion are kept, namely, sampling of high vorticity regions, of downflow zones, and reduction of the average rise velocity [9].
The comparison of cases with and without lift force and equal v T , τ b (i.e., cases (i) and (ii) in table 2) shows that the respective diffusion coefficients are basically the same. When τ b = τ k case (iii), the clustering in vortex structures is enhanced and the rise velocity is reduced: there is a competition between stronger trapping, that should reduce the diffusion, and slowing down, that enhances it. The second effect prevails and the diffusivities are larger than before. Note also that for τ b = τ k , the fluctuations of all observables are much stronger than in the other cases, because of the strong interaction of bubbles with the small scale turbulence structures.
The picture which arises is that bubbles are trapped within small scale eddies and move upwards together with them until they break up. At this time they are released. The trapping reduces the dispersion during the life-time of the structure. Once released the fast decay of the correlation between the fluid velocity at the position of the eddy break up and at the initial bubble position further prevents an efficient dispersion. We conclude that the interaction with vortices, due to the flow pressure gradient in the bubble motion equation (i.e., the 3Du/Dt term), is the main origin of the large discrepancy between numerical results and analytical theory.
The bubble diffusivity D b i is linked to the integral timescale T b,i L via the velocity fluctuations (v b i ) 2 , see equation (23). We computed in all simulations vertical and lateral bubble velocity fluctuations of the order of u 2 0 , independently of β. This indicates that the bubble integral timescale 2 follows the same behaviour of D b i , as function of β. The result agrees with analytical predictions for µ → ∞, which corresponds to negligible inertia effects, i.e., motion governed by drag and gravity.
In simulations with two-way coupling the difference between vertical and lateral dispersion is smaller. Indeed, we showed that the effect of bubbles on turbulence [9] is the attenuation of large-scale vertical fluid velocity fluctuations. Thus the diffusion is reduced in the gravity direction. At the same time, it may be enhanced in the horizontal planes, due to the increased velocity fluctuations in these directions (caused by incompressibility, see table 3 of [8]); however, this effect is small, as seen from teh nearly agreeing curves for the horizontal diffusivities in figure 2.
flow structures (that occurs here because τ b τ k ), bubble motion resembles the one of fluid particles.
We remark that, according to our results, the maximum diffusivity achieved is equal to the fluid particle diffusivity. In SB97, the authors claim the existence of a blow up of bubble velocity fluctuations at large β, and a subsequent unlimited increase of the diffusivity Our calculations do not support this statement. However, in that work large β values were attained by stronger turbulent fluctuations u 0 , while keeping v T and µ = L 11 /(τ b v T ) constant. These constraints induce an increase of the bubble response time τ b compared to the turbulence timescales. In such regime, viscous forces are weak and cannot counteract the rapidly varying inertial forces. The consequence might be the stated blow up of the fluctuations (v b i ) 2 . In our case we change v T , whereas the bubble timescale is fixed and is always smaller than all turbulence timescales (τ b = τ k /10), thus viscous forces are always dominant in determining the bubble velocities.
In figures 4 and 5, the lateral and vertical Lagrangian velocity autocorrelation functions R i L (t) are plotted for several values of β. The data refer to simulations with one-way coupling. The curves indicate that for smaller β, the lateral velocity decorrelates faster in time. For β = 0.7, the function has even a negative loop that explains the small diffusivity achieved in this regime (see figure 2). The crossing trajectory effect is a reasonable explanation for the observed trend in β. The fluid particle velocity autocorrelation is an upper limit for these curves.
On the other hand, there is no pronounced trend in β for the vertical autocorrelation, as long as β 1, i.e., as long as the bubble-rise velocity is smaller than the r.m.s. fluid velocity fluctuations, the turbulent structures can trap the bubbles [36]. In contrast, for β = 0.7, the autocorrelation drops to zero faster than before. The results suggest that it is the combination of bubble trapping by turbulent structures that occurs when v T < u 0 (i.e., β > 1) and of the preferential bubble collection in the downflow side of these structures that leads to a vertical bubble diffusivity of the same order of the fluid particle one. However, the latter autocorrelation is not an upper limit to R z L (t).

Relative dispersion
To further explore the differences between fluid particles and bubbles it is interesting to compare their relative diffusions that measure how the distance between two initially close particles increases as a function of time [2,37]. The relative dispersion tensor for either fluid particles or bubbles is defined according to: where Y i (t) is the ith component of the separation vector between two particles: with separation Y i0 at the initial time t = t 0 : We fix here t 0 = 0. The brackets · indicate average over all particle pairs. The presence of an average motion in the z-direction does not change the definition of the bubble dispersion tensor in comparison to the one of fluid particles, because we are concerned with relative distances between bodies that are all translating upwards. , v T = 3v k (· · · · · ·), and v T = 6v k (---). For comparison the fluid relative dispersion tensor is also presented (--). The straight dotted lines indicate the scaling in the viscous, in the Richardson, and in the time asymptotic regimes.
Three regimes for the relative dispersion tensor are to be expected in the large Re regime [38,39]: (i) at small separation time t the behaviour of the tensor is quadratic in time, more precisely with v rel,i the ith component of the relative velocity v rel = v 1 (t) − v 2 (t). (ii) For large enough flow Reynolds numbers, i.e., in the case of fully developed turbulence (which is not yet given in our simulations with Re λ = 62), a regime for intermediate timescales, called Richardson regime, exists in which the dispersion is proportional to the third power of time [2,39,40] with the energy flux.
Therefore (see equation (20)), it holds for the diagonal components of the relative dispersion tensor i.e., the diagonal tensor components are linear in time with a proportionality coefficient that is twice the one that applies to absolute dispersion.
In figures 6 and 7, the relative dispersions of bubbles in horizontal planes and in the gravity direction are presented. Indeed, all three regimes (i), (ii), and (iii) can be recognized. For comparison, the dispersion of fluid particles is also plotted, which is not very different.   The tensor is computed on pairs of bubbles/fluid particles whose separation at t = 0 is between 2η and 4η.
At small time separations, all curves show agreement with dispersion of fluid particles (solid lines in the figures). Indeed, bubbles are released with initial velocity equal to the one of the surrounding fluid, thus in the beginning they follow the fluid particle trajectories. Later the behaviour changes and gradually fits to the asymptotic regime, in which the relative intensity of dispersion depends on β, as we showed in the previous paragraph (see figure 2 and equation (30)).
In figures 8 and 9, the logarithmic local slopes of the relative dispersions are plotted as functions of time. The graphs indicate that the asymptotic scaling t 3 is never attained, however there is a clear trend from the power 2 scaling at small t to a steeper one at larger t and then again a decrease up to the linear scaling for large time intervals.

Pdfs of the forces
The motion of fluid particles and bubbles, and therefore their absolute and relative dispersion, is determined by forces that, in turbulence, display a high degree of intermittency. We present here the pdfs of these forces. Bubble forces, as well as fluid particle accelerations, are evaluated at all mesh points and then interpolated out of grid, at the bubble position, by employing the third-order Taylor series scheme (TS13) [32].

Pdf of forces acting on fluid particles
In figure 10, we present the pdf of the fluid acceleration evaluated along the trajectories of fluid particles. For comparison, the pdf(F A ) along the bubble trajectories is also shown. The figure indicates that for fluid particles the tails of the curve are slightly higher than for bubbles, thus the intermittency is more pronounced.

Pdf of forces acting on bubbles
We examine now the other forces acting on bubbles. The effect of drag and buoyancy is given by the force Here we substracted the buoyancy to allow for the comparison between the three directions. For the lift force it holds Results for the pdfs are presented in the simulation with one-way coupling (i.e., τ b = τ k /10 and v T = 2v k ). Figures 11-13 show that the pdfs are strongly intermittent, and this property holds especially for the lift force, computed in the z-direction, that probes a very intermittent quantity, namely the small scale vorticity field.
In table 3, the flatnesses f = F X 4 /( F X 2 ) 2 of the forces F X acting on bubbles are presented. For comparison, the flatness of the acceleration computed along fluid particle trajectories is also shown. Once again we remark the strong intermittency displayed by the lift force.
with F X the fluid acceleration, F A , the drag force, F D (defined in equation (31), and the lift force F L , evaluated along the bubble trajectories. For comparison, in the last column, the values measured for the fluid acceleration along fluid particle trajectories are also reported. Results are shown in the horizontal planes (xy) and in the vertical direction (z) separately. We may infer from the results that, if one considers a regime in which bubble motion is ruled by inertial forces rather than by viscous forces, sudden bursts in the intensity of the lift can lead to a strong increase of the bubble velocity fluctuations. This regime is the one studied by SB97. Therefore, such reasoning can justify the blow up of bubble velocity fluctuations that is employed in that paper to explain the divergence of the diffusivity at high β values.

Lagrangian velocity structure functions
Up to now in literature the main focus has been on Eulerian fluctuations and on the strong intermittency displayed by spatial velocity increments [43]. However, more recently, intermittency has been measured experimentally also for Lagrangian fluctuations [41]- [45] by tracking the temporal evolution of fluid particle velocities and accelerations. The numerical method employed here offers the possibility to investigate and to compare Lagrangian statistics for fluid particles and for bubbles.
As we showed in the previous section, the acceleration of both of them is a highly intermittent quantity. The same property holds for the velocity increments τ v i (t), evaluated at small time lags τ: where v i (t) is the ith component of the Lagrangian velocity of a particle or bubble and τ the time interval. The pdfs of τ v i (t) for fluid particles are shown, for several time intervals τ, in figure 14. For the smallest τ the pdf is as wide as the one of the acceleration (also reported in the figure), whereas at large τ it can be well approximated by a Gaussian function.

Lagrangian statistics for fluid particles
In the Lagrangian approach, the velocity structure function of order p is defined by The average is over time t. In the inertial range, when τ k τ T L , with τ k and T L , respectively, the Kolmogorov and the integral Lagrangian timescales, according to Kolmogorov's 1941 hypothesis, the second-order moment should only depend on the energy flux and on τ, Here C 0 is a universal constant. Recent experimental measurements [41] as well as numerical analysis [46]- [48] converge in the range C 0 = 6 ± 0.5. In figure 15 D L 2 (τ), normalized by the product τ is presented as a function of τ for the present numerical simulation. Owing to the low Reynolds number, no plateau is observed. However, the curve reaches the maximum at C * 0 ∼ 3, which is in agreement with the predicted trend at small Reynolds numbers [49], C * 0 /C 0 = 0.07Re 1/2 λ ≈ 0.55. For the higher order moments, the corresponding Kolmogorov dimensional prediction is with ζ(p) = p/2. As in the Eulerian approach to turbulence, any deviation of the scaling exponent ζ(p) from the linear law are related to intermittency. In the Lagrangian framework, the intermittency is even stronger than in the Eulerian, as shown within the multifractal approach and also through direct numerical simulations [50]. It can be quantified by comparing the measured values of the scaling exponents to the Kolmogorov predictions. However, a direct measure of ζ(p) is not possible due to the low Reynolds number. The inertial range in which the scaling law (35) holds is in general not wide enough. Therefore, just for the study of the Eulerian scaling  figure 18) from τ/T L = 0.25 to τ/T L = 2. Note that according to Kolmogorov's dimensional prediction ξ(p) = p/2 (· · · · · ·).
properties, it is common to employ the extended self similarity (ESS) approach [51], that consists of plotting the pth order structure function against the qth order structure function, p = q. Such plots normally show reasonable scaling, even for structure functions that do not display scaling when plotted versus the time interval τ, see e.g., choose figure 15. We q = 2 and obtained a scaling relation D L p (τ) = (D L 2 (τ)) ξ(p) , with ξ(p) = ζ(p)/ζ (2). The structure functions of order p = 1, 3, 4 and 5 resulting from our simulations are presented in figure 17, as functions of the second order one. The straight lines that fit the data have slopes, from top to the bottom, 0.56, 1.34, 1.56 and 1.8, in good agreement to the experimental results [41] shown in figure 16. The logarithmic local slopes of D L p (τ) versus D L 2 (τ) are plotted in figure 18.
In figure 16, we plot all the exponents ξ(p) = ζ(p)/ζ (2) for the Lagrangian moments, as obtained from our numerics. For comparison we also show the scaling exponents for the Eulerian structure functions, as obtained by the She-Leveque model [52], that closely fits to turbulence measurements at high Reynolds numbers. The Lagrangian intermittency corrections are much larger than the Eulerian ones. The experimental results [41] are also presented in the figure; they nicely agree with our numerical findings. The numerical findings for ξ(p) of the simulations presented in [50] are slightly larger, but also agree with the experimental results and our numerical results within the error bars.

Lagrangian statistics for bubbles
The Lagrangian velocity structure functions can be evaluated also along the bubble trajectories. Thus,  as in equation (33), but now the velocity v is the bubble's instantaneous velocity. Scaling properties and intermittency issues can be addressed also for these structure functions. We focus here on three cases: (i) v T = 2v k and τ b = τ k /10; (ii) v T = 4v k and τ b = τ k /10; and (iii) v T = 6v k and τ b = τ k /10. Only one-way interactions are considered.
In figure 19, the second-order structure functions are presented for case (iii). We distinguish between the Lagrangian statistics in the gravity (z) direction and in the horizontal (xy) planes. In fact, owing to the difference between the Lagrangian timescales, T b,xy L < T b,z L (see figure 2 and the relation between diffusivity and timescale (23)), the possible scaling ranges are not expected to coincide. The figure displays that, indeed, the xy structure function saturates much faster with time than the z structure function. For comparison, the fluid particle structure function is also plotted.
As for the fluid Lagrangian statistics, we assess scaling by using the ESS method. In figure 20, we present the logarithmic local slope of D L 5 (τ) versus D L 2 (τ), for xy-and z-directions separately, as well as for the total function. The plot indicates that it is difficult to detect any scaling exponent, because the inertial range is not well developed. However, it looks like that there is a transition between two different laws for the structure function in the z-direction. The scaling exponent increases from a value close to 1.8 to another one close to 2.2, which represents, respectively, the Lagrangian and the Eulerian intermittent scaling exponent. The transition occurs at τ/T L ∼ 0.8.
A possible recovery of the Eulerian scaling in the gravity direction can be explained by looking at the vertical velocity autocorrelation function. When v T u 0 and µ = L 11 /(τ b v T ) 1, the bubble velocity fluctuations can be approximated by the fluid velocity fluctuations evaluated along the bubble trajectory [5,53] with v T ∝ (0, 0, v T ). Now, the local fluid velocity fluctuations occur on a timescale ∼L 11 /u 0 , that is much larger (as β is small) than the timescale ∼L 11 /v T on which the fluctuations due to bubble motion occur. Thus one can drop the Eulerian time dependence in equation (37) which indicates that the Lagrangian velocity autocorrelation in the z-direction can be approximated with the two point spatial correlation of Eulerian velocities. With this argument in mind, it should no longer be too surprising that we detect similar scaling properties for bubble Lagrangian and fluid Eulerian statistics. In figure 21, the logarithmic local slope for structure functions along the gravity direction is presented for all three cases (i) to (iii), with 0.7 < β < 2. The graph depicts that, with the decrease of β, the local slope seems to increase from the Lagrangian to the Eulerian scaling law in the same way as the Lagrangian velocity autocorrelation function tends towards the Eulerian one.

Conclusions
The dispersion of microbubbles with response time τ b much smaller than the Kolmogorov time τ k has been studied in homogeneous isotropic turbulence and compared to the dispersion of fluid particles. The numerics has been compared with the analytical predictions of SB98, who address the same bubble regime, namely, µ = L 11 /(v T τ b ) 1 and β = u 0 /v T ∼ 1, with L 11 , v T and u 0 , the integral turbulence scale, the bubble-rise velocity in still liquid, and the fluid velocity r.m.s. fluctuations, respectively. Since τ b τ k bubble motion is mainly ruled by viscous forces. However, inertial forces are fundamental in the dispersion: the numerical results deviate from the analytical theory mainly owing to the interaction with eddies, which is caused by the flow pressure gradient term in the bubble equation of motion. In particular, the lateral diffusivities (in directions perpendicular to the gravity) are remarkably smaller than what is predicted in the theory. This work further supports the importance of bubble clustering in dispersion processes. Such analysis can only be achieved by a full direct simulation of the flow, i.e., including a proper representation of the small scales, which is not provided by hitherto performed kinematic simulation. The forces acting on bubbles are strongly intermittent, and this holds in particular for the lift force that probes a very intermittent quantity, namely, the small scale vorticity field. The analysis of Lagrangian intermittency is carried out by measuring the velocity structure function along both fluid particle and bubble trajectories. The numerical results confirm, as already claimed in the past [41], that intermittency is stronger in the Lagrangian framework than in the Eulerian. For bubbles, at large β values, we find similar intermittent scaling laws as for fluid particles. However, at small β, the presence of a strong drift along gravity makes the bubble Lagrangian velocity autocorrelation in this direction close to the spatial correlation of fluid Eulerian velocities. In this case, also the intermittency displayed by the bubble velocity structure functions is closer to the fluid Eulerian intermittency rather than to the fluid Lagrangian.
The effect of an artificial (homogeneous and isotropic) forcing that substains the turbulence reflects in the numerical values of the diffusion coefficients. However, relative values of bubble-to-fluid particle diffusivities, as well as the qualitative trends of dispersion tensors and pdfs, are linked to turbulence characteristics and thus to the inertial range of scales. As a consequence, the results are not expected to qualitatively depend on the form of the external forcing that drives and affects only the large flow scales. The situation would be different for anisotropic or inhomogeneous forcing. No predictions on those cases can be drawn from our work.