Preferential concentration vs. clustering in inertial particle transport by random velocity fields Piero Olla ISAC-CNR and INFN, Sez. Cagliari, I–09042 Monserrato, Italy. (Dated: January 13, 2010) Theconceptofpreferentialconcentrationinthetransportofinertialparticlesbyrandomvelocity fields is extended to account for the possibility of zero correlation time and compressibility of the 0 velocity field. It is shown that,in thecase of an uncorrelated in time random velocity field,prefer- 1 entialconcentrationtakestheformofaconditiononthefieldhistoryleadingtothecurrentparticle 0 positions. This generalized form of preferential concentration appears to be a necessary condition 2 for clustering in the uncorrelated in time case. The standard interpretation of preferential concen- n tration is recovered considering local time averages of the velocity field. In the compressible case, a preferentialconcentrationoccursinregionsofnegativedivergenceofthefield. Intheincompressible J case, it occurs in regions of simultaneously high strain and negative field skewness. 3 1 PACSnumbers: 05.10.Gg,05.40.-a,46.65.+g ] h I. INTRODUCTION approach, in the case of inertial particle transport, has c allowedderivationofanalyticalexpressionsforthe parti- e m Oneofthemoststrikingcharacteristicsofinertialpar- cle concentration correlations, both in the weak [17, 20] ticle transport by random velocity fields is clustering. and large [21] inertia limits. As recognized in [17], the - t This phenomenon occurs, but is not confined to, in tur- weakinertialimit,inaKraichnanmodelapproach,corre- a t bulent flows [1]; clustering phenomena, in fact, were ini- sponds to a regime of adiabatic variationfor the particle s tially predicted to occur with particles pushed by a 1D separation. In this framework, the particle concentra- . t (one-dimensional)randomforcefield[2]. Theinteresting tion dynamics can be cast in the form of a problem of a m pointisthataninitialspatiallyhomogeneousdistribution fast variable elimination, in which the fast variables are the particle velocities. of inertial particles will develop clumps and voids, even - d iftheflowisincompressible. Bothexperimentalevidence A mechanism that has been proposed for cluster for- n [1] and numerical simulations [3, 4] confirm this effect. mation in turbulence, is the centrifugal force induced, o Spatialinhomogeneity ofthe randomfield statistics may preferentialconcentrationofheavy(light)particlesinthe c contribute, but is not crucial to the process. strain (vorticity) regions of the fluid [22]. This preferen- [ It is to be mentioned that clustering phenomena are tialconcentrationeffectwaslaterconfirmedinnumerical 4 thought to be important both in industrial flows [5], in simulations [3, 23]. Now, clustering turns out to occur v the atmosphere (the problem of rain formation) [6, 7] alsoin1D andinthe Kraichnanmodel just discussed,in 7 and in the oceans (the problem of plankton dynamics, situations therefore, where it is not clear what meaning 2 especially as regards blooming) [8, 9]. shouldbegiventopreferentialconcentration. Inparticu- 9 2 Overtheyears,asubstantialtheoreticalefforthasbeen lar,theconceptsofstrainandvorticitydonotexistin1D. . directed to the analysis of clustering by random flows This casts some doubts on whether preferential concen- 1 [10–17]. As first noticed in [2], clustering appears to be tration (or some generalized version of it) is an essential 1 8 associated with a weak inertia regime, in which the sep- ingredientforinertialparticleclusteringinrandomflows. 0 aration of the trajectories of the particles and the fluid In general, preferentialconcentration could be defined : elements they crossin their motion, is small onthe scale as the fact that averages of fluid (random flow) quanti- v of the trajectory evolution. This corresponds to a situa- ties, obtained from sampling along inertial particles tra- i X tion in which the relaxation time of the particle velocity jectories, do not coincide with what would be obtained r isshorterthanthecharacteristictimeoftherandomfield from spatial (or temporal) averages. In other words, it a fluctuations. In the absence of molecular diffusion, the could be interpreted as a non-ergodicity property of the resulting clusters are ofsingular nature,concentratedon process of random field sampling by the particle flow. a set of zero measure [2, 18]. Given the situation, a first question that would be in- Anapproachthathasbeenfruitfulinthestudyofpas- terestingto ask,is whether the clustering ofinertialpar- sive tracer transport, is that of considering uncorrelated ticles in random flows is always the result of the non in time random velocity fields, the so called Kraichnan ergodic sampling, by the particles, of some relevant field model [19]. The role of characteristic time of the ran- quantity. Inthepresenceoffinitecorrelationtimeandin- dom field is played in this case by the diffusion time of compressibility, physical considerations allowed to iden- a tracer (or, depending on the problem, pairs of tracers) tify from the start, strain and vorticity as the relevant across a correlationlength of the field. Weak and strong quantities. In the general case, such an operation may inertia will then refer to fast or slow relaxation with re- be not as easy, and an interesting question is, therefore, spect to this characteristic time. The Kraichnan model whether the correct relevant quantities could be identi- 2 fied directly from the equations of motion. to frozen field conditions. For real turbulence, K ∼ 1 Thispaperwilltrytoanswerthesetwoquestions. The and S parameterizes the strength of inertia. analysiswillshowthattheanswerispartoftheprocedure The K → 0 regime, corresponding to a Kraichnan offastvariableeliminationrequiredtodeterminethedy- model for inertial particles, allows to neglect the par- namics of the particle concentration field. It is actually ticle displacement in a correlation time τ . Clustering E the way in which the fast variable elimination procedure is a phenomenon associated with singular behavior, at handles the presence of memory in the original process. small values of the argument, of the particle separation Thispaperisorganizedasfollows. InSec. II,theequa- PDF (probability density function) ρ(r,t). A Kraichnan tions for the transport of a pair of inertial particle in a model approach allows to describe the two-particle dy- randomfieldarederived. InSec. III,evolutionequations namics relevant for clustering, in terms of a single equa- for the particle separation distribution are derived and tion: applied to the determination of the clustering strength. ν˙ =uˆ(r,t)−ν, r˙ =ν, (5) In Sec. IV the issue of preferential concentration and its contribution to clustering is discussed. Section V is astheevolutionofuˆ(r,t)dependsonlyonr,andnotsep- devoted to conclusions. Technical details are left in the aratelyonthecoordinatesofthetwoparticles. [Forfinite Appendices. τ , this would not be true, as uˆ(r(t),t) would depend, E on scale τ , on the separate evolution of the two parti- E cle coordinates, described by Eq. (1)]. In a Kraichnan II. MODEL EQUATIONS modelapproach,theseparationPDFρ(r,t)willobeythe Fokker-Planck equation associated with Eq. (5) (more Themotionofaninertialparticleinarandomvelocity precisely, its restriction to the variable r). field u(x,t) can be modelled by the Stokes equation As recognized in [16], the two-particle dynamics be- comes dependent, in an uncorrelated in time regime, on v˙ =τ−1[u(x,t)−v], x˙ =v. (1) the single parameter S Influidmechanics,thismodelwoulddescribethedynam- ǫ= σu2τE =K2S, (6) ics of a particle that is sufficiently small, of sufficiently r2 v high density [24], and in flow conditions in which the ef- that is the amplitude factor in front of Eq. (2). This fect of gravitycan be neglected. The relaxationtime τ , S parameter plays the role of generalized Stokes number calledthe Stokestime, inthe caseofasphericalparticle, for a Kraichnan model. Writing ǫ = τ /τ˜ , one notices isgivenbyτ =2/9a2λ/ν ,withaistheparticleradius, S E S 0 in fact that τ˜ = r2/(σ2τ ) is the time for a tracer [a λ the ratio of the particle and fluid densities and ν the E v u E 0 point moving with velocity u(x(t),t)] to diffuse across fluid kinematic viscosity. r , that plays the role of characteristictime scale for the Consider a smooth, D-dimensional Gaussian veloc- v uncorrelated in time random flow. ity field u(x,t), with stationary and spatially uni- Smallǫcorrespondstoaweakinertiaregime,andEqs. form and isotropic statistics. For an incompressible (2,5,6) provide, for τ → 0, a Kraichnan model for the random field, the structure function huˆ (r,t)uˆ (r,0)i, E α β clusteringofweaklyinertialparticles(see[16,17,20]and uˆ(r,t)≡u(x+r,t)−u(x,t), can be written in the form references therein). To see that small ǫ corresponds to huˆ (r,t)uˆ (r,0)i=F(t)gˆ (r), where, for r ≪r : α β αβ v weak inertia, it suffices to verify that r changes little on a time τ (the correlation time for ν): ∆r(τ ) ≪ r. σ2τ S S gˆ (r)= u E (D+1)r2δ −2r r , (2) This is verified a-posteriori solving Eqs. (5,2) for fixed αβ rv2 h αβ α βi r, r ≪ r . This gives hν2|ri ∼ ǫr2, that coincides, in v ∞ the dimensionless units of Eq. (4), with the square dis- and 0 |F(t)|dt = 1, with σu, rv and τE, respectively, placementinatime τ . Thus,∆r(τ )/r ∼ǫ−1/2 ≪1,as the cRharacteristic velocity, length and time scales of the S S expected. field. Following [16], the following dimensionless param- This is an adiabatic regime for r, that has allowedthe eters are introduced: authors in [17] to use a fast variable elimination tech- nique to derive a version of the Fokker-Planck equation S =τS/τE, K =σuτE/rv, (3) associated with Eq. (5), restricted to r. In the following sections, a similar approach will be utilized to establish called respectively Stokes and Kubo numbers. Units are the connection between clustering and preferential con- chosen such that σ =τ =1; in this way: u S centration effects. τE =S−1, rv =(KS)−1. (4) The Kubo number describes the intrinsic long- or short- III. CLUSTERING correlated in time nature of the field, with K → 0 (S → ∞) corresponding to an uncorrelated regime: Thefastvariableeliminationprocedureinastochastic F(t) → 2δ(t); K → ∞ (S → 0) corresponds in turn problem, like the one described by Eq. (5) in the ǫ ≪ 1 3 regime, can be carried on substantially in two ways [25]: Following the approach in [29], outlined in Appendix working at the Fokker-Planck equation level, by means A, the response function can be calculated as an expan- of so called projection operator techniques [26, 27], or sion in powers of uˆ [i.e., from Eqs. (2) and (6), basically averaging away the fast variables already at the level of inpowersofǫ1/2]: R (t;z,τ)=δ(r(τ)−z)[Rˆ(0)+Rˆ(1)+ γβ γβ γβ the stochastic differential equation. The procedure fol- ...]. The first two terms in the expansion read [see Eqs. lowed in [17] was of the projection operator type, along (A3-A4)]: the lines of the approach derived in [28], in the context ofstochasticclimatemodelling. Thisapproachisnotap- Rˆ(0) = ψ(τ −t)δ , ψ(t)=θ(−t)[1−et], (10) γβ γβ propriate here: dealing with preferential concentration t ∂uˆ (r(τ′),τ′) effectswillrequireevaluationofconditionalaveragessuch Rˆ(1) = dτ′ψ(τ −τ′)ψ(τ′−t) γ (11) as hfn[uˆ]|r(t)=¯ri, with fn[uˆ] products of the fields and γβ Zτ ∂rβ(τ′) their derivatives, about which, the Fokker-Planck equa- andθ(t)istheHeavisidestepfunction[θ(t)=1fort>0, tion associated with Eq. (5) provides no information. θ(t)=0 otherwise]. The approach that is going to be followed here, is the From Eqs. (10-11) and (A4), it appears that onedescribedin[29,30],basedontheuseofthesocalled R (τ;z′,τ) = 0, so that the functional derivative on stochasticLiouvilleequation,andoffunctionalderivation γβ δ(r(τ) − z), arising in Eq. (7) from uˆ (r(τ),τ) = techniques [31, 32]. A similar approach was followed in α dDzuˆ (z,τ)δ(r(τ)−z), doesnotcontribute. Theanal- [21], to analyze the large inertia limit of particle trans- α Rysis in Appendix A shows that the lowest order inertia port. contributiontoJ occursatO(ǫ2),correspondingtotak- The evolutionequationfor the PDF ρ(¯r,t), associated ing into accountαonly R(1)(t;z′,τ) in the expansion for with Eq. (5), can be written in the form ∂ ρ(¯r,t) + γβ ∂¯ J (¯r,t) = 0, ∂¯ ≡ ∂/∂r¯ , with the probabtility cur- Rγβ(t;z′,τ). Carrying out the necessary functional dif- α α α α ferentiationsinEq. (7)andusingEqs. (10-11)givesthen rent J given by α the result, for t=0: t J (¯r,t)= dτeτ−thuˆ (r(τ),τ)δ(r(t)−¯r)i. (7) 0 α Z α J ≃ −2∂¯ dτeτ gˆ (r(τ))δ(r(0)−¯r) δ ψ(τ) −∞ α γZ−∞ D αβ h γβ This evolution equation is basically the average over all 0 ∂uˆ (r(τ′),τ′) ′ ′ ′ γ tfihgeurraeatiloiznatsipoancseooffuˆE,qo.ft(h5e):L∂ioρu˜(v¯ri,llte)e+qu∂¯at[iνon(tin)ρ˜t(h¯re,tc)o]n=- + Zτ dτ ψ(τ −τ )ψ(τ ) ∂rβ(τ′) iE. (12) t α α 0,ρ˜(¯r,t)=δ(r(t)−¯r). Theintegral t dτeτ−tuˆ(r(τ),τ) Thelowestordercontributiontothecurrentis,exploiting inEq. (7),infact,isjustthesolutioRn−∞forν(t)ofEq. (5), incompressibility, ∂¯γgˆαγ(¯r)=0: and J≡hν(t)|r(t)=¯riρ(¯r,t), with ρ(¯r,t)=hρ˜(¯r,t)i. J(2) =−g (¯r)∂¯ ρ(¯r,t), (13) The correlation between the Dirac delta and the ran- α αγ γ dom field in Eq. (7) is calculated using the functional corresponding,inEq. (12),to putgˆ (r(τ))≃gˆ (r(0)) αβ αβ integration by part formula [31, 32] (see also [33]). Indi- andtoneglecttheterm∝∂ uˆ intheintegral. Equation β γ cating by δ/δuˆ(x,t) the operation of functional differen- (13) describes pure tracer transport. tiation, this corresponds to making in Eq. (7) the sub- The next order is J(4), that receives two contribu- stitution tions: the one ∝∂ uˆ , originating from Rˆ(1), and one ∝ β γ gˆ (r(τ))−gˆ (r(0)). Thefirstcontributioniscalculated δ αβ αβ uˆ (z,t)−→2 dDr′gˆ (z,z′) , (8) applying Eq. (8) to the factor ∂uˆ (r(τ′),τ′)/∂r (τ′) α Z αβ δuˆ (z′,t) γ β β in (12). The second contribution is calculated writing where 2gˆ (z,z′)δ(t − t′) = huˆ (z,t)uˆ (z′,t′)i and of gˆαβ(r(τ))−gˆαβ(r(0))asafunctionof∆r(τ)=r(τ)−r(0) course gˆαβ(z,z) ≡ gˆ (z) [34].α In orβder to use Eq. and using the expression, which descends from Eq. (5): αβ αβ (8), however, it is first necessary to write in Eq. (7) 0 uˆ (r(τ),τ) = dDzuˆ (z,τ)δ(r(τ) − z); the functional ∆r(τ)= dτ′Φ(τ′,τ)uˆ(r(τ′),τ′), (14) α α Z derivative δ/δuˆR (z′,τ), will then act on a product −∞ β δ(r(τ) − z)δ(r(t) − ¯r). This requires determination of where expressions like ′ ′ ′ Φ(τ ,τ)=ψ(τ −τ)−ψ(τ ). (15) δr (t) Rγβ(t;z′,τ)= δuˆ (γz′,τ). (9) EvaluatingtoO(ǫ2)thetwocontributionsandexploiting β incompressibility, leads to the result (see Appendix B): Writing δ/δuˆ (z′,τ)δ(r(0)−¯r) = −R (t;z′,τ)∂¯ δ(r(t) β γβ γ 1 −¯r), it is clear that the role of the response function J(4) =− ρ(¯r)∂¯ ∂¯ gˆ (¯r)∂¯ gˆ (¯r) (16) R (t;z′,τ) in Eq. (7) is precisely to account for the α 2 β η αγ γ βη γβ correlation between the random field and the condition plus other terms containing spatial derivatives of ρ. The on the separation at time t. stationary solution ρ(¯r) is obtained requiring that the 4 current be divergence free: ∂ J = 0. The lowest order creasing powers of uˆ. The first contributions are: α α solution, from Eq. (13), is spatially uniform, and the derivatives acting on ρ do not contribute in Eq. (16), at f1[uˆ] = ∂¯αuˆα(¯r,τ), (21) the order considered. The stationary solution will obey, f [uˆ] = ∂¯ ∂¯ [uˆ (¯r,τ)uˆ (¯r,τ′)], (22) 2 α β α β therefore, the equation f [uˆ] = ∂¯ ∂¯ ∂¯ [uˆ (¯r,τ)uˆ (¯r,τ′)uˆ (¯r,τ′′)]. (23) 3 α β γ α β γ gˆ (¯r)∂¯ ∂¯ ρ(¯r)+ ρ(¯r)∂¯ ∂¯ gˆ (¯r)∂¯ ∂¯ gˆ (¯r)=0, (17) Other contributions involve gradients of the Dirac delta αγ α γ β η αγ α γ βη 2 in Eq. (7), that would lead in the end to gradients of ρ(¯r,t). The terms that would lead to clustering even thatisinagreementwith[17],onceawrongsignintheir starting from a spatially homogeneous particle distribu- Eq. (40) is corrected [35]. As discussed in [17, 20], the tion, however, are those in Eqs. (21-23). solution of Eq. (17), with the expression for gˆ (r) pro- αβ The contribution from Eq. (21) is evaluated following vided in Eq. (2), is a power law at r → 0: ρ(r) ∝ r−cǫ, the same procedure leading from Eq. (7) to (12). Stop- c = 2(D +1)(D +2), corresponding to the correlation ping to lowest order in ǫ and neglecting terms involving dimensionfor the particle distribution: D =D−2(D+ 2 gradients of ρ, leads to the result (see Appendix C): 1)(D+2)ǫ. 1 h∂¯ uˆ (¯r,τ)|r(0)=¯ri=− ψ(τ)∂¯ ∂¯ gˆ (¯r), (24) α α α γ αγ 2 IV. PREFERENTIAL CONCENTRATION that will vanish if u is incompressible. This means that, in the case of a compressible random field, particularly In the previous section, clustering was obtained by inD =1,clusteringwillbetheresultofpreferentialcon- solving an evolution equation for the separation PDF centration in regions of negative divergence of u. More ρ(¯r,t): ∂tρ +∂¯αJα = 0, in which the probability cur- precisely, particle pairs that at time t = 0 have separa- rent divergence ∂¯αJα was expressed, through Eq. (7), in tion ¯r, in the past were more likely to be in regions of termsofaconditionalaverage∂¯αhuˆα(r(τ),τ)δ(r(0)−¯r)i. negative ∇·u. Since, in the absence of conditioning, this averagewould ItistobenoticedthatEq. (18)doesnotcorrespondto be zero, it is clear that non-ergodic sampling of the ran- anexpansionofJ in powersof ǫ1/2. The fields uˆ in Eq. α dom field along the particle trajectories is a necessary (18), upon substitution in Eq. (7), contribute O(ǫ1/2) in condition for a non-zero current, and hence for cluster- contraction with another field uˆ, but contribute O(ǫ) in ing. contractionwithδ(r(0)−¯r). Aconsequenceofthisisthat In order to clearly identify the quantities undergoing Eq. (24) does not provide the whole O(ǫ) part of ∂¯ J . α α preferential concentration, it is necessary to pass the di- In the compressible case, another O(ǫ) contribution to vergenceoperatorin∂¯αhuˆα(r(τ),τ)δ(r(0)−¯r)i,insidethe ∂¯αJα isprovidedfor instance by∂¯α∂¯βgˆαβ(¯r), thatcomes average,sothattheexpressioncanbewrittenintheform from the part of h∂¯ ∂¯ [uˆ (¯r,τ)uˆ (¯r,τ′)]|r(0)=¯ri that is α β α β hf[uˆ]δ(r(0)−¯r)i. From f[uˆ], one can then extract the uncorrelated with δ(r(0)−¯r). contributions from specific monomials fn[uˆ] in the fields In the incompressible case, the contribution from f1 and their derivatives. to the current divergence vanishes, and it is necessary In the present situation, a good strategy to identify to consider f , with n > 1; in particular, the term in n thesemonomialsistoexpandthefielduˆ (r(τ),τ)around Eq. (22). Now, the strain andthe vorticityvariance of a α the final point r(0): velocity field u are defined from uˆ(r(τ),τ) ={1+∆rβ(1)(τ)∂β +∆r(2)(τ)∂β |S|2 =(1/2)∂αuβ(∂αuβ+∂βuα) 1 |ω|2 =∂ u (∂ u −∂ u ), (25) + ∆r(1)(τ)∆r(1)(τ)∂ ∂ +...}uˆ(r(0),τ), (18) α β α β β α 2 β γ β γ and ∂ u ∂ u = |S|2 − |ω|2/2, whose average can be α β β α where shown to be zero if u is incompressible. Thus, f [uˆ] is 2 connected with the difference |S|2−|ω|2/2 for u. ∆r(1)(τ) = Φˆuˆ(r(0),τ), (19) Since, in the case of an incompressible field: h|S|2i = ∆r(2)(τ) = Φˆ∆r(1)∂ uˆ(r(0),τ), (20) h|ω|2/2i,hf2[uˆ]|r(0)=¯riwillreceiveitsfirstnon-zerocon- β β tributionbycontractionofthetwofieldsuˆ anduˆ inf α β 2 with δ(r(0)−¯r) [see Eq. (22)]. The calculation carried andsoontohigherorders,withΦˆg(τ)≡ 0 dτ′Φ(τ′,τ) ′ −∞ on in Appendix C gives the result, for τ >τ : ×g(τ′) and Φ(τ,τ′) given in Eq. (15). InRthe weak iner- tia regime considered, ǫ ≪ 1, this is appropriate in the h∂¯ uˆ (¯r,τ′)∂¯ uˆ (¯r,τ)|r(0)=¯ri=ψ(τ) α β β α interval −τ ≤1, in which the time integral in Eq. (7) is ×[3ψ(τ′)−ψ(τ′−τ)]∂¯ ∂¯ gˆ (¯r)∂¯ ∂¯ gˆ (¯r), (26) α β γη γ η αβ concentrated. Substituting Eqs. (18-20) into (7), the current diver- plus terms involving gradients of ρ(¯r,t). The quantity gence can be written as a sum of averages, involving in- to RHS (right hand side) of Eq. (26) is positive for 5 ′ τ <τ <0. Thus,preferentialconcentrationinregionsof erential concentration, through an expansion in the ef- highstrain,inthewidersensegivenheretotheterm,isa fective Stokes number ǫ [see Eqs. (3,4,6)]. In the two property that is mantained in an incompressible Kraich- regimesof compressibleand incompressible flow,the rel- nan model for inertial particle transport. evant field properties are related to the field divergence Again, Eq. (26) does not account for the whole of the at the particle position, in the first case, and to the dif- current divergence to O(ǫ2) in the incompressible case. ference between the square strain and vorticity in the One has to consider also the skewness-like term in Eq. second. An additional relevant quantity in the incom- (23). This time, it is necessary to consider the contrac- pressiblecase,notpreviouslyconsideredinthe literature tion of only one of the fields with δ(r(0)−¯r). A calcula- (to the author’sknowledge),is the skewnessin Eq. (27). tion completely analogousto that of Eq. (24) gives then The standard interpretation of preferential concentra- the result: tion is recovered considering a local time average of the flow: u¯(¯r,t) = ∆t−1 ∆t/2 uˆ(¯r,t+τ)dτ. Basically, this h∂¯ ∂¯ ∂¯ [uˆ (¯r,τ)uˆ (¯r,τ′)uˆ (¯r,τ′′)]|r(0)=¯ri −∆t/2 α β γ α β γ gives to the randomRfield a finite correlation time ∆t. =−∂¯α∂¯βgˆγη(¯r)∂¯γ∂¯ηgˆαβ(¯r)[ψ(τ)δ(τ′−τ′′) From Eqs. (21-23), it is then possible to define time av- +permutations], (27) eraged quantities ∂¯αu¯α(¯r,τ), ∂¯α∂¯β[u¯α(¯r,τ)u¯β(¯r,τ)], and ∂¯ ∂¯ ∂¯ [u¯ (¯r,τ)u¯ (¯r,τ)u¯ (¯r,τ)], corresponding to diver- α β γ α β γ plus, again, terms involving gradients of ρ(¯r,t). Thus, gence, difference between strain and vorticity squared clustering receives, in the incompressible case, con- [seeEq. (25)]andskewnessofthefieldu¯. FromEqs. (24- tribution from preferential concentration in regions 27), it is easy to see that inertial particles, with respect of simultaneously high strain and negative skewness tothesequantities,undergopreferentialconcentrationin ∂¯ ∂¯ ∂¯ [uˆ (¯r,τ)uˆ (¯r,τ′)uˆ (¯r,τ′′)]. the usual sense of the word. α β γ α β γ The terms in Eqs. (26-27) are all the preferential con- The analogy between particle transport by a Kraich- centration contributions that arise, at O(ǫ2), in the in- nan model and by a random velocity field with finite compressible case. A term f of fourth order in uˆ, not correlation time, with the parameter ǫ in the first case 4 indicated in Eq. (18), may contribute to hf [uˆ]|r(0)=¯ri; set equal to S in the second, allows predictions of the 4 at O(ǫ2), this occurs through the unconditional average preferential concentration strength in the finite correla- hf [uˆ]i. However,Gaussianityofureducesthiscontribu- tion time case. Indicating by h.i , average at a parti- 4 p tion to products hf2[uˆ]ihf2[uˆ]i, which vanish in the case cle position: (rv/σu)h∇ · uip = O(S1/2), in the com- of an incompressible random field. pressible case, and (r /σ )2h(|S|2 −|ω|2/2)i = O(S), v u p It is to be noticed that the RHS of Eqs. (24-23) are (r /σ )3h∂ ∂ ∂ [u u u ]i = O(S1/2) in the incompre- independent of the argument¯r. This means that the av- sibvle uone. αPβreγfereαntβialγcopncentration should occur in eragesinthoseequationsareconditionedtothepresence regions of negative field divergence in the compressible of a pair of particles at unspecified separation. In other case, high strain and negative skewness in the incom- words,the two-particleconditionalaveragesin Eqs. (24- pressible one. 27), are equivalent to averages conditioned to a single The techniques utilized in this paper are not specific particle in x+¯r. to inertial particles, and are in fact rather general. For instance, conclusions analogous to the ones on inertial particles, could be drawn, in the case of fluid tracers in V. CONCLUSION compressibleflows,bothregardingclusteringandprefer- ential concentration. (In fact, it would be interesting to The analysis in this paper shows that some general- understandhow much of particle clustering in compress- ized form of preferential concentration continues to be ibleflowsisassociatedwithinertia,andhowmuchisdue an ingredient for clustering, also in situations where its to accumulation on fluid shocks). presence is not expected, such as transport by uncorre- lated in time, or 1D random velocity fields. The original statement, that the random field at the particle location Acknowledgments has on the averagesome property (e.g. higher strain), is replaced, in the uncorrelated case, by one on properties TheauthorwishestothankK.Turitsynforinteresting of the random field, still at that location, but at times and helpful discussion. previous to the arrival of the particle. The current field configuration, instead, is uncorrelated with the current particle positions. In other words, if the field is uncorre- Appendix A. Response function determination latedintime, therewillbe nopreferentialconcentration, in the standard meaning of the term, rather, a general- Following [29], an equation for the response function izedversioninvolvingfieldandparticle configurationsat R (t;z,τ) = δr (t)/δuˆ (z,τ) can be obtained writing γβ γ β different times. Eq. (5) in the form The weak inertia considered, has allowed to identify the quantities involved in this form of generalized pref- r¨ (t)+r˙ (t)=uˆ (r(t),t) γ γ γ 6 and taking the functional derivative with respect to plus terms involving spatial derivatives of ρ. uˆ (z,τ). The result is, for t>τ: The other contribution in Eq. (12) is: β R¨γβ(t;z,τ)+R˙γβ(t;z,τ)= ∂uˆ∂γ(rrη((tt)),t)Rηβ(t;z,τ). Jαg = −2∂¯γZ 0 dτeτψ(τ) −∞ (A1) FromEq. (5),onecanwriter(t)=const.+ t dτψ(τ− × h[gˆαγ(r(τ))−gˆαβ(r(0))]δ(r(0)−¯r)i, (B3) −∞ t)uˆ(r(τ),τ) and ν(t) = t dτeτ−tuˆ(r(τ),τR), with ψ(t) which can be seen to contain both quadratic and linear −∞ given in Eq. (10). ThesRe expressions lead to the initial terms in ∆r, with ∆r given in Eq. (14). Writing conditions gˆ (r(τ))−gˆ (r(0))=Gˆ [∆r (τ)∆r (τ) R (τ;z,τ)=0; R˙ (τ;z,τ)=δ δ(r(τ)−z). (A2) αβ αβ αβηφ η φ γβ γβ γβ +r (0)∆r (τ)+r (0)∆r (τ)], η φ φ η Equation (A1) with the initial conditions in Eq. (A2) can be solved perturbatively in uˆ (i.e. basically in ǫ1/2): with Gˆαβηφ = (1/2)∂η∂φgˆαβ(r), substituting into Eq. R (t;z,τ)=δ(r(τ)−z)[Rˆ(0)+Rˆ(1)+...],withtheresult (B3) and using Eq. (14) allows to calculate immediately γβ γβ γβ the quadratic term Rˆ(0) = ψ(τ −t)δ , (A3) γβ γβ 5 Rˆ(n+1) = tdτ′ψ(τ′−t)∂uˆγ(r(τ′),τ′)Rˆ(n). (A4) Jαg,2 =−6ρ(¯r)∂¯β∂¯ηgˆαγ(¯r)∂¯γgˆβη(¯r), (B4) γβ Z ∂r (τ′) ηβ τ η plus terms involving, again, spatial derivatives of ρ; the In particular, the second order term reads: linear terms lead only to contributions involving deriva- tives of ρ. Combining Eqs. (B2) and (B4) leads to Eq. t τ′ Rˆ(2) = dτ′ dτ′′ψ(τ −τ′′)ψ(τ′′−τ′)ψ(τ′−t) (16). γβ Z Z τ τ ∂uˆ (r(τ′),τ′)∂uˆ (r(τ′′),τ′′) γ φ × ∂r (τ′) ∂r (τ′′) . (A5) Appendix C. Evaluation of non-ergodic terms φ β It is possible to see that, to O(ǫ2), this term does not As with Eq. (7), calculation of conditional averages contribute to the current Jα. In fact, including Rˆγ(2β) in like hfn[uˆ]|r(0)=¯ri, with fn[uˆ] given in Eqs. (21-23), is Eq. (12), the contraction∂ uˆ ∂ uˆ ∝ǫδ(τ′−τ′′) would carriedonbyapplicationofthefunctionalintegrationby φ γ β φ be killed by the factor ψ(τ′′ −τ′) from Eq. (A5). On part formula Eq. (8). In the case of f [uˆ], one has to 1 the other hand, the contractions of the fields ∂ uˆ and lowest order in ǫ, from Eqs. (8) and (10): φ γ ∂ uˆ , with the other factors in Eq. (12), would pro- dβuceφadditionalfactorsǫ1/2,beyondthosefromthefields h∂¯αuˆα(¯r,τ)δ(r(0)−¯r)i themselves, and the final result would be of higher order = 2[∂¯ hδ(r(0)−¯r)∂¯ gˆ (¯r,r(τ))i β α αβ in ǫ. − hδ(r(0)−¯r)∂¯ ∂¯ gˆ (¯r,r(τ))i] (C1) α β αβ From the definition: 2gˆ (r,r′)δ(t − t′) = αβ Appendix B. Contributions to the probability huˆ (r,t)uˆ (r′,t′)i, with uˆ(r,t) = u(x+r,t)−u(x,t), it α β current is easy to see that gˆ (r,r′) = (1/2)[gˆ (r)+gˆ (r′)− αβ αβ αβ gˆ (r−r′)], and therefore αβ The contribution from Rˆ(1) in Eq. (12) is evaluated using the functionalintegrationbypartformulaEq. (8), ∂rαgˆαβ(r,r′)|r′=r =(1/2)∂rαgˆαβ(r). (C2) and keeping only Rˆ(0) in the resulting expansion for the SubstitutingintoEq. (C1)withtheapproximation(valid response function. Indicating by JR this contribution: α to the order considered) gˆ (¯r,r(τ))≃gˆ (¯r), leads im- αβ αβ 0 0 mediatelytoEq. (24). Thecalculationofhf [uˆ]|r(0)=¯ri, JR = 4 dτ dτ′ dDz dDz′eτψ(τ −τ′) 3 α Z Z Z Z that leads to Eq. (27), is completely analogous. −∞ τ The calculation of hf [uˆ]|r(0)=¯ri is slightly more in- × ψ2(τ′)∂ gˆ (z,z′)∂¯ ∂¯ hgˆ (r(τ)) 2 zβ γη γ η αβ volved and requires considering the correction terms in × δ(r(τ′)−z)δ(r(τ′)−z′)δ(r(0)−¯r)i. (B1) the response function, Eq. (11). The starting point is the equation ′ Notice that, due to the condition τ > τ, the functional derivativeactedonlyonδ(r(0)−¯r)andnotongˆαβ(r(τ)). huˆ (¯r,τ)uˆ (¯r,τ′)δ(r(0)−¯r)i=4ψ(τ)ψ(τ′) To O(ǫ2), one sets r(τ′)=r(0) in Eq. (B1) and carrying ×αhgˆ (¯r,βr(τ)gˆ (¯r,r(τ′))∂¯ ∂¯ δ(r(0)−¯r)i out the integrals obtains the result αγ βη γ η −2huˆ (r(τ),τ)gˆ (¯r,r(τ′))R(1)(0;r(τ′),τ′) 1 α βγ γφ JαR = 3ρ(¯r)∂¯β∂¯ηgˆαγ(¯r)∂¯γgˆβη(¯r), (B2) ×∂¯φδ(r(0)−¯r)i. (C3) 7 Forτ >τ′,theproductuˆ (r(τ),τ)R(1)(0,r(τ′),τ′)inthe τ)ψ(τ)∂¯ gˆ (¯r,r(τ)). Substituting into Eq. (C3), ap- α ηφ γ αη second term to RHS of Eq. (C3), from Eq. (11), leads proximating gˆαβ(¯r,r(τ)) ≃ gˆαβ(¯r,r(τ′)) ≃ gˆαβ(¯r) and to a factor 2ψ(τ′ −τ)ψ(τ)∂ gˆ (¯r,r(τ)) = 2ψ(τ′ − using Eq. (C2), leads, after little algebra, to Eq. (26). rγ(τ) αη [1] J.R.Fessler,J.D.KulickandJ.K.Eaton,Phys.Fluids6, [18] J. Bec, Phys. Fluids 15, L81-L84 (2003) 3742 (1994) [19] R.H. Kraichnan, Phys.Rev.Lett. 72, 1016 (1994) [2] J.M. Deutsch,J. Phys. A:Math. Gen. 18, 1449 (1985) [20] M. Wilkinson, B. Mehlig, S. Ostlund and K.P. Duncan, [3] R.C.HoganandJ.N.CuzziPhys.Fluids13,2938(2001) Phys. Fluids 19, 113303 (2007) [4] G. Boffetta, F. De Lillo and A.Gamba Phys. Fluids 16, [21] P. Olla and M.R. Vuolo, Eur. Phys.J. B 65, 279 (2008) L20 (2004) [22] M.R. Maxey,J. Fluid Mech. 174, 441 (1987) [5] C.T. Crowe, M. Sommerfeld and Y. Tsuji Multiphase [23] K.-P. Wang and M.R. Maxey, J. Fluid Mech. 256, 27 Flows with Particles and Droplets (CRC Press, New (1993) York,1998), and references therein. [24] M.R. Maxey and J.J. Riley,Phys. Fluids 26, 883 (1983) [6] P.A.VaillancourtandM.K.Yau,Bull.Am.Met.Soc.81, [25] K.LindenbergandB.J. West The nonequilibrium statis- 285 (2000) tical mechanics of open and closed systems (VCH Pub- [7] R.A.Shaw, Annu.Rev. Fluid Mech. 35, 183 (2003) lishers, 1990) [8] J. Ruiz,D. Macias and F. Peters, Proc. Natl. Acad. Sci. [26] R. Zwanzig, Phys. Rev.124, 983 (1961) U.S.A.101, 17720 (2004) [27] P. Grigolini and F. Marchesoni, Adv. Chem. Phys. 62, [9] R. Reigada, R.M. Hillary, M.A. Bees, J.M. Sancho and 29 (1985) F. Sagues, Proc. R. Soc. Lond.B 270, 875 (2003) [28] A.Majda, I.Timofeyev and E. VandenEijnden, Comm. [10] T. Elperin, N. Kleeorin and I. Rogachevskii, Phys. Rev. Pure Appl.Math. 54, 891 (2001) Lett.77 5373 (1996) [29] J.M.Sancho,M.SanMiguel,S.L.KatzandJ.D.Gunton, [11] E. Balkovsky, G. Falkovich and A. Fouxon, Phys. Rev. Phys. Rev.A 26, 1589 (1982) Lett. 86 2790 (2001) [30] P. H¨anggi, T.J. Mroczkowski, F. Moss and P.V.E. Mc- [12] T. Nishikawa, Z. Toroczkai and C. Grebogi, Phys. Rev. Clintock, Phys.Rev. A 32, 695 (1985) Lett.87, 038301 (2001) [31] K. Furutsu, J. Res. Nat. Bur. Stand. Sect. D 67, 303 [13] S.Sigurgeirsson and A.M.Stuart,Phys.Fluids14, 4352 (1963) (2002) [32] E.A. NovikovSov.Phys. JETP 20, 1290 (1965) [14] J.C.H.FungandJ.C.Vassilicos,Phys.Rev.E68,046309 [33] U. Frisch, Turbulence: the legacy of A.N. Kolmogorov (2003) (Cambridge University Press, 1995) [15] L.I.ZaichikandV.M.Alipchenkov,Phys.Fluids15,1776 [34] For a given realization, the two velocity differences (2003) uˆ(r,t)=u(x+r,t)−u(x,t)anduˆ(r′,t′)=u(x+r′,t′)− [16] K. Duncan, B. Mehlig, S. O¨stlund and M. Wilkinson, u(x,t′) are defined at the same point x (the choice of x Phys.Rev.Lett. 95, 240602 (2005) is part of therealization selection). [17] J. Bec, M. Cencini, R. Hillerbrand and K. Turitsyn, [35] K. Turitsyn,private communication Physica D 237, 2037 (2008)