Heavy particles in a persistent random flow with traps

TL;DR

This paper models heavy particles in a one-dimensional compressible fluid with persistent Gaussian velocity fields, analyzing trapping, Lyapunov exponents, and caustic formation, with implications for higher-dimensional systems.

Contribution

It introduces a novel analytical framework for understanding particle trapping and caustic formation in a 1D compressible flow with weak inertia.

Findings

Particles tend to accumulate in trapping regions.

Analytical expressions for Lyapunov exponents and caustic rates are derived.

Numerical simulations confirm the analytical predictions.

Abstract

We study a one-dimensional model for heavy particles in a compressible fluid. The fluid-velocity field is modelled by a persistent Gaussian random function, and the particles are assumed to be weakly inertial. Since one-dimensional fluid-velocity fields are always compressible, the model exhibits spatial trapping regions where particles tend to accumulate. We determine the statistics of fluid-velocity gradients in the vicinity of these traps and show how this allows to determine the spatial Lyapunov exponent and the rate of caustic formation. We compare our analytical results with numerical simulations of the model and explore the limits of validity of the theory. Finally, we discuss implications for higher-dimensional systems.

Figures

Equations

$$\dot{x}(t)=v(t),\dot{v}(t)={\tau }_{\text{p}}^{-1}\left[-v(t)+u(x(t),t)\right],$$

$$\displaystyle\langle u(x,t)u(x^{\prime},t^{\prime})\rangle=u_{0}^{2}\big{[}1\!-\!(x\!-\!x^{\prime})^{2}/\eta^{2}\big{]}\text{e}^{-(x-x^{\prime})^{2}/\eta-|t-t^{\prime}|/\tau_{\text{c}}}\,.$$

$$\text{Ku}={\tau }_{\text{c}}{\tau }_{\text{a}}^{-1},\text{St}={\tau }_{\text{p}}{\tau }_{\text{c}}^{-1}\text{and}\kappa =\sqrt{3}{\tau }_{\text{p}}{\tau }_{\text{a}}^{-1}.$$

$$\dot{X}$$

$$\dot{\bar{x}}$$

$$\dot{Z}={\tau }_{\text{p}}^{-1}\left[A(\bar{x},t)-Z\right]-Z^2.$$

$$\dot{X}$$

$$\dot{Z}$$

$$P(A_n=a)=\{\begin{array}{ll}\frac{|a|}{{\kappa }^2}{\text{e}}^{-a^2/(2{\kappa }^2)}& a\le 0,\\ 0& \text{otherwise}.\end{array}$$

$$⟨A_n⟩$$

$$⟨A_n{A}_m⟩-⟨A_n⟩⟨A_m⟩$$

$${\mathcal{F}}_{A_n=a}\varrho (Z=z|A_n=a)={\partial }_t\varrho (Z=z|A_n=a),$$

$${\mathcal{F}}_{A_n=a}=\frac{\text{d}}{\text{d}z}(z+z^2-a).$$

$$\frac{\text{d}}{\text{d}z}(z+z^2-a){\varrho }_s(Z=z|A_n=a)=0.$$

$${\varrho }_s(Z=z|A_n=a)=\{\begin{array}{ll}\frac{\sqrt{-1-4a}}{2\pi \left(z^2+z-a\right)}& a<-1/4\\ \delta (z-z^{\ast })& a\ge -1/4\end{array},$$

$$P_s(Z=z)=$$

$$P_s$$

$$+{\int }_0^{\infty }\frac{\sqrt{a}(a+1/4)}{\pi {\kappa }^2[(z+1/2{)}^2+a]}{\text{e}}^{-(a+1/4{)}^2/(2{\kappa }^2)}\text{d}a,$$

$$\lambda =\underset{t\to \infty }{\lim }\underset{X(0)\to 0}{\lim }\frac{1}{t}\log \frac{|X(t)|}{|X(0)|},$$

$$\lambda$$

$$\lambda =-\frac{1}{2}+\frac{{}_2{F}_2\left(1,\frac{3}{2};\frac{7}{4},\frac{9}{4};-\frac{1}{32{\kappa }^2}\right)}{120{\kappa }^2},$$

$$\mathcal{J}=\underset{z\to \infty }{\lim }{z}^2{P}_s(Z=z)={\int }_0^{\infty }\text{d}a\frac{\sqrt{a}(a+1/4)}{\pi {\kappa }^2}{\text{e}}^{-(a+1/4{)}^2/(2{\kappa }^2)}.$$

$$\mathcal{J}=(2^{\frac{5}{2}}\pi {)}^{-1}{\text{e}}^{-1/(64{\kappa }^2)}{K}_{-\frac{1}{4}}\left(\frac{1}{64{\kappa }^2}\right),$$

$$\lambda =\underset{t\to \infty }{\lim }\frac{1}{t}{\int }_0^{t}Z(s)\text{d}s.$$

$$\mathcal{J}=\underset{t\to \infty }{\lim }\frac{N(t)}{t}.$$

$$\ddot{X}+\dot{X}-A_{n}X=0.$$

$$X(t)=C_1{\text{e}}^{{\lambda }_{+}t}+C_2{\text{e}}^{{\lambda }_{-}t},$$

$$\lambda$$

$$=-\frac{1}{2}+{\int }_{-1/4}^{\infty }\text{d}a\sqrt{1/4+a}P(A_n=a).$$

$$\mathcal{J}$$

$$=\frac{1}{\pi }{\int }_{-\infty }^{-1/4}\text{d}a\sqrt{-1/4-a}P(A_n=a),$$

Full text

I Introduction

Turbulent fluids in Nature and Technology often contain small particles whose density is larger than that of the carrier fluid. Examples are sprays, water droplets in turbulent clouds, and dust in the exhaust of combustion engines. It is a common observation that these particles do not distribute homogeneously in space, but instead form regions of high and low concentration on small scales. This effect is referred to as spatial clustering Bec (2003); Wilkinson and Mehlig (2003); Gustavsson and Mehlig (2016), and it is a consequence of the interaction of the particles with the turbulent fluid. Hence, any attempt to describe these systems must take into account the turbulent fluctuations, a challenging problem by itself.

Recently, much progress has been made in explaining and quantifying spatial clustering of heavy particles in turbulence, using statistical models. See Ref. Gustavsson and Mehlig (2016) for a recent review. Such models approximate the turbulent fluid-velocity fluctuations below the smallest scales of turbulence, called the dissipative range, in terms of random velocity fields with prescribed statistics. Clearly, these random velocity fields lack the complexity of real turbulence. Remarkably, however, such statistical models often reproduce the dynamics of the immersed particles and the statistical properties of their spatial distribution qualitatively, and in some cases even quantitatively. Statistical models allow for much simpler numerical simulations and even admit analytical treatment in limiting cases Wilkinson and Mehlig (2003); Gustavsson and Mehlig (2016); Falkovich et al. (2007). The models typically depend on two parameters, the Kubo number Ku, which is a measure of the persistence of the flow, and the Stokes number St, which measures the importance of particle inertia in the fluid-particle system.

Spatial clustering is often characterised in terms of the spatial Lyapunov exponents $\lambda_{i}$, $i=1,\ldots,d$ Sommerer and Ott (1993); Gustavsson and Mehlig (2016). These exponents measure the rates at which particles approach or diverge, in the long-time limit. Because heavy particles are not constrained to follow the flow, their dynamics takes place in $2d$-dimensional phase space, while their spatial concentration is measured in the $d$-dimensional coordinate space. As a consequence, the spatial distribution of the particles exhibits singularities, so-called caustics, which contribute to the spatial clustering of the particles Falkovich et al. (2002); Wilkinson and Mehlig (2005); Wilkinson et al. (2006); Falkovich et al. (2007). One therefore identifies the rate of formation of these singularities, the rate of caustic formation $\mathcal{J}$, as a second important observable.

A well-studied limiting case of the statistical model is the so-called white-noise limit, describing the dynamics of very heavy particles in a rapidly fluctuating fluid-velocity field. In terms of the dimensionless parameters, this limit corresponds to $\text{Ku}\to 0$ and $\text{St}\to\infty$ so that $\text{Ku}^{2}\text{St}$ remains constant Gustavsson and Mehlig (2016). An important property of white-noise models is that they depend solely on a single parameter, $\varepsilon^{2}=f(d,\beta)\,\text{Ku}^{2}\text{St}$, where $f$ is a function of the spatial dimension $d$ and the compressibility parameter $\beta$ of the underlying flow field. The parameter $\varepsilon$ takes the role of an inertia parameter in the white-noise model, similar to the Stokes number St in statistical models at finite values of Ku Gustavsson and Mehlig (2016). In the white-noise limit in one spatial dimension, $\lambda=\lambda_{1}$ and $\mathcal{J}$ can be computed analytically as a function of $\varepsilon$ Wilkinson and Mehlig (2003).

In this Paper, we describe a second limiting case of the statistical model, the persistent-flow limit. This limit describes weakly inertial, heavy particles in a very persistent flow, hence the name. The limit is taken by letting $\text{Ku}\to\infty$ and $\text{St}\to 0$ such that the product Ku St stays finite. Similarly to the white-noise limit, the persistent-flow limit turns out to depend on a single inertia parameter $\kappa$, where $\kappa^{2}=(d+2)\text{Ku}^{2}\text{St}^{2}$ Gustavsson and Mehlig (2016).

Here, we confine ourselves to the one-dimensional version of the model, $d=1$, and show that $\lambda$ and $\mathcal{J}$ can be computed analytically in the persistent-flow limit. In one spatial dimension, the random fluid-velocity field is compressible. Because of the high persistence of the flow, compressibility leads to long-lived trapping regions, regions in space where the particles accumulate. We show that the statistics of fluid-velocity gradients in the trapping regions are the crucial ingredient to solving the one-dimensional model.

The Paper is organised as follows: In Sec. II we introduce the model for heavy particles in turbulence and motivate the persistent-flow limit by a detailed analysis of the relevant time scales. Section III contains the statement and the derivation of the analytical results for the spatial Lyapunov exponent $\lambda$ and the rate of caustic formation $\mathcal{J}$. We compare these results to those obtained in the white-noise limit, and to results of numerical simulations of the statistical model. We draw conclusions and discuss implications for higher-dimensional models in Sec. IV. The Appendix contains an alternative derivation of the main formulae derived in Sec. III.

II Problem formulation

We study a one-dimensional model for the motion of heavy particles in the dissipative range of turbulence, approximating the viscous force on the particle by Stokes drag. The equations of motion for an underdamped small and heavy particle in a viscous fluid are given by

[TABLE]

where $\tau_{\text{p}}=\gamma^{-1}$ is the particle relaxation time, the typical time it takes the particle to relax to the fluid velocity. We model the turbulent fluid velocity $u(x,t)$ by a one-dimensional Gaussian random velocity field with zero mean and correlation Gustavsson and Mehlig (2016)

[TABLE]

In this formulation, the flow field is a function of its root mean square velocity $u_{0}$ as well as its correlation length and fluid correlation time, $\eta$ and $\tau_{\text{c}}$, respectively. Fields in one spatial dimension are always potential, therefore the velocity field $u(x,t)$ is compressible, as mentioned in the Introduction. Models for particles in three-dimensional turbulence, by contrast, usually assume that the flow is incompressible. Further, approximating $u(x,t)$ by a Gaussian random field neglects several characteristic aspects of turbulence such as intermittency and broken time-reversal invariance. It has, however, been shown that the highly-simplified one-dimensional model gives important insights into the dynamics of particles in the dissipative range of turbulence (see Gustavsson and Mehlig (2016) and references therein).

II.1 Relevant time scales

The random velocity field $u(x,t)$ is characterised by two different time scales. In addition to the correlation time $\tau_{\text{c}}$, the advection time scale $\tau_{\text{a}}=\eta/u_{0}$ is the typical time it takes an advected (tracer) particle to travel one correlation length $\eta$. From the three time scales $\tau_{\text{c}}$, $\tau_{\text{a}}$ and $\tau_{\text{p}}$, one can define three dimensionless numbers Gustavsson and Mehlig (2016)

[TABLE]

The Kubo number, Ku, measures the persistence of the flow and St is the Stokes number, a measure for particle inertia in the system. The third quantity in expression (3), $\kappa$, compares the particle relaxation time to the advection time and is thus another measure of particle inertia.

In case there is a time-scale separation between $\tau_{\text{a}}$ and $\tau_{\text{c}}$, only one of the two inertia parameters, either St or $\kappa$, needs to be considered. In the white-noise limit Wilkinson and Mehlig (2003, 2005); Wilkinson et al. (2006); Derevyanko et al. (2007); Gustavsson and Mehlig (2016); Pumir and Wilkinson (2016), for instance, the correlation time $\tau_{\text{c}}$ is the smallest of all time scales. More precisely, this limit is approached if $\tau_{\text{p}}\gg\tau_{\text{a}}\gg\tau_{\text{c}}$ so that the parameter $\varepsilon\propto\sqrt{\tau_{\text{p}}\tau_{\text{c}}}\tau_{\text{a}}^{-1}={\rm Ku}\sqrt{{\rm St}}$ stays finite. In other words, the white-noise limit of Eqs. (1) describes very heavy particles in a rapidly fluctuating flow field. The fluid-velocity gradient becomes a white-noise signal with zero mean and correlation strength $\propto\varepsilon^{2}$. Thus, Eqs. (1) constitute a generalised diffusion process and one can use Fokker-Planck equations to solve the model exactly in one dimension, $d=1$ Wilkinson and Mehlig (2003); Gustavsson and Mehlig (2016). The single relevant dimensionless inertial parameter in the white-noise model is $\varepsilon$, which is related to the diffusion constant $D$ of the generalised diffusion by $D\tau_{\text{p}}^{-1}\propto\varepsilon^{2}\propto\text{Ku}^{2}\text{St}$.

In the persistent-flow limit studied here, by contrast, the correlation time $\tau_{\text{c}}$ is the largest of the three time scales, i.e., $\tau_{\text{c}}\gg\tau_{\text{p}},\tau_{\text{a}}$. In terms of dimensionless parameters this limit corresponds to $\text{Ku}\gg 1$ and $\text{St}\ll 1$. The ratio between the remaining time scales is thus given by $\kappa\propto\tau_{\text{p}}/\tau_{\text{a}}$. Hence, the particle motion and the motion of advected particles take place on time scales much smaller than the correlation time of the fluid. The inertia parameter $\kappa$ measures how much the particles detach from the flow, compared to the motion of advected tracers.f

The advection time scale $\tau_{\text{a}}$ in the statistical model is analogous to the Kolmogorov time $\tau_{\text{K}}$ in homogeneous isotropic turbulence, defined as $\tau_{\rm K}=1/\sqrt{{\rm tr}\,\langle\mathbb{A}\mathbb{A}^{\sf T}\rangle}$. Here $\mathbb{A}$ is the matrix of fluid-velocity gradients, and $\langle\cdots\rangle$ is an average along fluid-element paths. Since the Kolmogorov time is the only available microscopic time scale in turbulence, the analogue of $\kappa$ (constructed from $\tau_{\text{K}}$ and $\tau_{\text{p}}$), is commonly used as the inertia parameter in DNS and experiments on the dynamics of heavy particles in turbulence Pergolizzi (2012); Gustavsson and Mehlig (2016).

II.2 One-dimensional model

The observables $\lambda$ and $\mathcal{J}$ are most conveniently studied in terms of the relative dynamics of an ensemble of particle pairs, instead of the single-particle dynamics given in Eqs. (1). The motion of a pair of particles is described by the particle separation $X=x^{(1)}-x^{(2)}$ and relative velocity $V=v^{(1)}-v^{(2)}$, as well as their midpoint position $\bar{x}=(x^{(1)}+x^{(2)})/2$ and mean velocity $\bar{v}=(v^{(1)}+v^{(2)})/2$. At separation much smaller than $\eta$, $X\ll\eta$, a particle pair experiences an essentially linear velocity field. We therefore linearise the equations of motion (1) with respect to $X$. To this end, we introduce the mean flow $\bar{u}(\bar{x},t)\sim\left[u(x^{(1)},t)+u(x^{(2)},t)\right]/2$ and the fluid-velocity gradient, $A(\bar{x},t)\sim\left[u(x^{(1)},t)-u(x^{(2)},t)\right]/X$. A particle pair at small separation then obeys the linearised equations of motion

[TABLE]

where we have omitted the time dependence of $(X,V,\bar{x},\bar{v})$ for convenience of notation. The particle-velocity gradient $Z=V/X$, has the equation of motion Falkovich et al. (2002); Wilkinson and Mehlig (2003); Gustavsson and Mehlig (2016)

[TABLE]

The two equations for $(\bar{x},\bar{v})$, Eqs. (4b), resemble a closed system that is equivalent to the single-particle dynamics (1). Equations (4a) and (5), on the other hand, are not closed as they depend on the variable $\bar{x}$ which enters in the argument of the fluid-velocity gradient.

As both advected tracers and heavy particles move much faster than the fluid changes, $\tau_{\text{a}},\tau_{\text{p}}\ll\tau_{\text{c}}$, and the fluid-velocity field $u(x,t)$ is compressible, there exist long-lived, stable fixed points of Eqs. (1), and thus Eqs. (4b), where the particles accumulate. These trapping regions are characterised by negative fluid-velocity gradients $A(x,t)<0$, and small fluid velocity $u(x,t)\ll u_{0}$. Trapping regions appear and disappear randomly with mean rate proportional to $\tau_{\text{c}}^{-1}$ Mehlig et al. (2009). When a trapping region vanishes, the trapped particles are driven to a neighbouring trapping region, at a typical distance $\sim\eta$ away.

If $\tau_{\text{p}}\ll\tau_{\text{a}}$ we estimate the migration time to the next trapping region by the time it takes for advected tracers to travel one correlation length, $\tau_{\text{a}}$. If $\tau_{\text{p}}\gg\tau_{\text{a}}$, on the other hand, the migration time is of the order of the relaxation time $\tau_{\text{p}}$ of the particles. We conclude that the migration time is of the order of the maximum of these two time scales, $\max(\tau_{\text{p}},\tau_{\text{a}})$.

As trapping regions persist for times $\sim\tau_{\text{c}}$, the persistent-flow limit $\tau_{\text{c}}\gg\tau_{\text{a}},\tau_{\text{p}}$ ensures that the time particles spend in trapping regions is much larger than the migration time. Hence, for $\text{Ku}\gg 1$ and $\text{St}\ll 1$ we can neglect the disappearance and appearance of stable zeros of $u(x,t)$ and consider the case of an ensemble of particles in an unbounded and frozen Gaussian random velocity field.

Each trapping region in the frozen velocity field is characterised by a constant fluid-velocity gradient $A_{n}<0$, where we index the trapping regions by $n$, for reasons that become clear later. Because different trapping regions are typically a fluid correlation length $\eta$ away from each other, their fluid-velocity gradients can be assumed to be independent and identically distributed.

To obtain the equation of motion for a particle pair in the $n^{\text{th}}$-trapping region of the flow, we first dedimensionalise Eqs. (4) and (5), so that they do not explicitly dependent on the time scales. That is, we let $t\to\tau_{\text{p}}t$, $X\to\eta X$, $V\to\eta\tau_{\text{p}}^{-1}V$, $A\to\tau_{\text{p}}^{-1}A$ and $Z\to\tau_{\text{p}}^{-1}Z$. We obtain the dimensionless equations

[TABLE]

In the rest of this Paper, all quantities, including the Lyapunov exponent $\lambda$ and the rate of caustic formation $\mathcal{J}$, are written in dimensionless form.

The probability density of the random fluid-velocity gradient in the $n^{\text{th}}$-trapping region, $A_{n}$, is given by the distribution of gradients at the zeros of the Gaussian random function $u(x,t)$, conditioned on that these gradients are negative Gustavsson and Mehlig (2013a). This distribution is obtained by means of the ‘Kac-Rice formula’ Kac (1948); Rice (1944), which allows to compute the density of singular points of random functions. Using that the gradient function $A(x,t)$ is Gaussian distributed, we obtain for the probability density of $A_{n}$:

[TABLE]

Eq. (7) shows, that the particle-velocity distribution depends only on the single parameter $\kappa$ [Eq. (3)] in the persistent limit. Furthermore it follows from Eq. (7) and the independence of the gradients $A_{n}$ that

[TABLE]

for all $n$ and $m$. Clearly, $\langle A_{n}\rangle<0$ is the result of strong preferential sampling of negative gradients, due to the accumulation of particles in trapping regions. Eq. (7) shows, in fact, that $P(A_{n}>0)=0$.

III Results

Using Eqs. (6) and the distribution of $A_{n}$ we calculate the spatial Lyapunov exponent $\lambda$ and the rate of caustic formation $\mathcal{J}$. Both quantities can be either obtained directly from the equations of motion for $X$ and $V$, Eqs. (6a), or from the steady-state probability density of $Z$, $P_{s}(Z=z)$. We explain the method using $P_{s}(Z=z)$ in the following Section and describe the first method in Appendix A. The general idea in both cases is to first compute the respective quantities conditional on a certain gradient realisation, $A_{n}=a$, and then average over all realisations using the probability density $P(A_{n}=a)$ given in Eq. (7).

III.1 Steady-state distribution of $Z$

We first turn to the computation of the steady-state distribution of $Z$. The conditional density $\varrho(Z=z|A_{n}=a)$ obeys the equation

[TABLE]

where the operator $\mathcal{F}_{A_{n}=a}$ is the generator of the Perron-Frobenius operator Cvitanovic et al. (2005) corresponding to the dynamics of $Z$ in Eq. (6b) given the realisation $A_{n}=a$. The operator reads

[TABLE]

By definition, the conditional steady-state probability density $\varrho_{s}(Z=z|A_{n}=a)$ solves $\partial_{t}\varrho_{s}(Z=z|A_{n}=a)=0$ and is thus annihilated by $\mathcal{F}_{A_{n}=a}$,

[TABLE]

Equation (11) has the solution

[TABLE]

where $z^{*}=1/2(-1+\sqrt{1+4a})$ is the stable fixed point of the dynamics (5) for $a>-1/4$. The drastic change of $\varrho_{s}(Z=z|A_{n}=a)$ at $a=-1/4$ can be understood by noting that Eq. (6b) undergoes a saddle-node bifurcation at that point, generating two fixed points of $Z$ for $A_{n}>-1/4$. For long times, $Z$ ends up at the stable of the two fixed points, $z^{*}$, in this case. For $A_{n}<-1/4$, on the other hand, there are no fixed points and $Z$ reaches $-\infty$ in a finite time. Because $Z=V/X$, this case corresponds to two particles overtaking each other, forming a caustic Gustavsson and Mehlig (2016); Falkovich et al. (2002); Wilkinson et al. (2006); Wilkinson and Mehlig (2005). Upon reaching $-\infty$, $Z$ is immediately reinjected at $+\infty$ to ensure $Z=V/X$ also after the caustic event Gustavsson and Mehlig (2011, 2014). We now calculate the steady-state probability density $P_{s}(Z=z)$, using

[TABLE]

and substituting the expressions in Eqs. (7) and (12) for $P(A_{n}=a)$ and $\varrho_{s}(Z=z|A_{n}=a)$, respectively. We obtain

[TABLE]

where, $\mathbf{1}_{S}$ denotes the indicator function of the set $S$. The steady-state distribution is the sum of two parts, corresponding to the two different solutions of $\varrho_{s}(Z=z|A_{n}=a)$. The first term in Eq. (14) originates from the gradient range $A_{n}\geq-1/4$, that does not admit caustics, while the second one comes from the gradient range where caustics form. The result for $P_{s}(Z=z)$, Eq. (14), as a function of $z$ for $\kappa=1/4$ is shown as the solid line in Fig. 1. The distribution exhibits a peak at $z=-1/2$ and a sharp kink at $z=0$. For large $z$ the distribution has power-law tails $\propto\mathcal{J}/z^{2}$, shown as the dotted lines. The steady-state distribution for $z$ in the persistent-flow model is thus quite similar to the corresponding quantity in the white-noise model, shown as the dash-dotted line in Fig. 1. We observe a similar peak at small $z$ and power-law tails with $z^{-2}$ scaling.

The apparent similarity of the two distributions in these very different limits is due to the similar qualitative dynamics of $Z$. In the persistent-flow model, $Z$ equilibrates at $z^{*}$, close to $-1/2$, for the majority of realisations, $A_{n}\geq-1/4$, hence the peak at $z=-1/2$. For $A_{n}<-1/4$, $Z$ repeatedly approaches $-\infty$ and returns, leading to the power law tails $\sim z^{-2}$ in the distribution. In the white-noise model, on the other hand, the noisy dynamics of $Z$ spends most of its time close to the stable fixed point at $z=0$. Occasionally, however, the gradient history allows $Z$ to pass the unstable fixed point at $z=-1$ and escape to infinity. The $z^{-2}$-tails of the distribution originate from caustic events $Z\to-\infty$. The large-$Z$ dynamics is universal and independent of the gradient $A(\bar{x},t)$, because the deterministic $Z^{2}$-term in Eq. (5) dominates for large enough $Z$.

III.2 Lyapunov exponent

The spatial Lyapunov exponent

[TABLE]

determines the long-time fate of particle separations. It is given by the steady-state expectation value $\lambda=\langle Z\rangle_{s}$ Gustavsson and Mehlig (2016):

[TABLE]

where $\mathcal{P}$ denotes the principle value of the integral. The second equality follows using Eq. (14) and performing integrations by part. The integral on the right-hand side of Eq. (16) is always negative and can be expressed in terms of special functions as

[TABLE]

where ${}_{2}F_{2}$ denotes the generalised hypergeometric function. The solid line in Fig. 2(a) shows the spatial Lyapunov exponent as a function of $\kappa$. For small $\kappa$, $\lambda$ scales as $\lambda\sim\langle A_{n}\rangle\propto-\kappa$, signalling clustering due to preferential concentration. For large $\kappa$, the spatial Lyapunov exponent asymptotes to $\lambda\sim-1/2+1/(120\kappa^{2})$ and thus approaches $-1/2$. The fact that $\lambda$ is negative for all $\kappa$ indicates that the particles in the trapping regions move closer and closer towards each other and the particle dynamics eventually contracts to a point for all $\kappa$.

As a comparison, the Lyapunov exponent for the white-noise model is shown as the solid line in Fig. 2(b). In the white-noise model, $\lambda$ is negative for small $\varepsilon$ and then changes sign and becomes positive at $\varepsilon_{\text{c}}\approx 1.33$. This behaviour is known as ‘path-coalescence transition’ Wilkinson and Mehlig (2003) or ‘aggregation-disorder transition’ Deutsch (1985) and describes the transition to chaotic behaviour of the model for $\varepsilon>\varepsilon_{\text{c}}$. For small $\varepsilon$, $\lambda$ scales as $\lambda\sim-\varepsilon^{2}$, meaning that $\lambda\propto\text{St}$ for small inertial parameter in both models. The scaling as a function of Ku is, however, different, as we observe $\lambda\propto\text{Ku}$ in the persistent-flow model and $\lambda\propto\text{Ku}^{2}$ in the white-noise model at small inertia parameter.

III.3 Rate of caustic formation

The rate of caustic formation $\mathcal{J}$ is the rate at which $|X(t)|\to 0$ while $|V(t)|>0$. That is, it is the rate at which particle trajectories cross at finite relative velocity. Globally speaking, it is the rate at which the phase-space manifold of the particles folds over. The rate of caustic formation is of particular importance for the calculation of collision rates between particles, where it leads to corrections to the Saffman-Turner formula for advected particles Saffman and Turner (1956); Pumir and Wilkinson (2016); Falkovich and Pumir (2007); Voßkuhle et al. (2014). We compute $\mathcal{J}$ from the magnitude of the tail of the steady-state distribution $P_{s}(Z=z)$, which behaves as $P_{s}(Z=z)\sim\mathcal{J}/z^{2}$ for large $z$:

[TABLE]

Again, the integral can be expressed in closed form using special functions. We obtain

[TABLE]

where $K_{\nu}$ is the modified Bessel function of the second kind. In Fig. 2(a) $\mathcal{J}$ is shown as the dash-dotted line. For small $\kappa$, $\mathcal{J}$ exhibits an exponentially small activation $\propto\text{e}^{-1/(32\kappa^{2})}$. Exponential activations are common for the rate of caustic formation of inertial particles in turbulence. They have been observed both in statistical models Gustavsson and Mehlig (2016); Wilkinson and Mehlig (2005); Wilkinson et al. (2006) and in numerical simulations Falkovich and Pumir (2007); Voßkuhle et al. (2014). However, the exponent of the activation differs, depending on the model and on the limit of Ku and St that is taken Gustavsson and Mehlig (2016); Mehlig and Wilkinson (2004); Gustavsson and Mehlig (2013b); Meibohm et al. (2017).

That said, the corresponding result for the white-noise model is shown as the dash-dotted line in Fig. 2(b). For small $\varepsilon$, $\mathcal{J}$ exhibits an activated behaviour $\propto\text{e}^{-1/(6\varepsilon^{2})}$. In terms of Ku and St, the activation has the exponents $-1/(96\,\text{Ku}^{2}\text{St}^{2})$ and $-1/(18\,\text{Ku}^{2}\text{St})$ for the persistent-flow model and the white-noise model, respectively.

III.4 Numerical simulations and deviations from persistent-flow limit

In order to study the validity of our analytical results for finite Ku and St, we perform direct numerical simulations of the statistical model. To this end, we generate a one-dimensional Gaussian random field with correlation (2) and let a large number of particles evolve according to the dynamics (1). The simulations are performed at large Ku ($10^{1},10^{2}$ and $10^{3}$), and different values of St, resulting in a range of values of the inertia parameter $\kappa=\sqrt{3}\,\text{Ku}\text{St}$.

The spatial Lyapunov exponent is obtained from direct numerical simulations by computing the long-time empirical mean of $Z$ along a particle trajectory $x(t)$ Gustavsson and Mehlig (2016):

[TABLE]

The rate of caustic formation is, in turn, calculated by counting the number $N(t)$ of the events $Z\to-\infty$ after a long time $t$, divided by $t$:

[TABLE]

Fig. 3(a) shows $\lambda$ obtained from Eq. (20) for $\text{Ku}=10^{1}$, $10^{2}$ and $10^{3}$ as pentagons, half circles and filled circles respectively. The solid line shows the analytic result (17). For small $\kappa$, the results of direct numerical simulations and the theory agree well for all measured Ku. This is expected, because for $\kappa\ll 1$ and fixed large Ku, the Stokes number is small, $\text{St}\ll 1$, and we are close to the persistent-flow limit. As Ku increases, the regime of validity for $\kappa$ increases and the theoretical result is approached. The persistent-flow theory says that $\lambda<0$ for all values of $\kappa$. We observe, however, that this prediction fails for finite (but large) values of Ku. Instead, $\lambda$ obtained from the direct numerical simulations using Eq. (20) exhibits a path-coalescence transition, just as in the white-noise model, for some $\kappa_{\text{c}}$. Hence, the dynamics becomes chaotic and the particles cluster on a fractal set for $\kappa>\kappa_{\text{c}}$, instead of contracting to a point.

Fig. 3(b) shows $\mathcal{J}$ obtained from direct numerical simulations of Eq. (21) for $\text{Ku}=10^{1}$, $10^{2}$ and $10^{3}$ as pentagons, half-circles and filled circles, respectively. The analytic result (19) is, as before, shown as the solid line. We observe similar behaviour as for the spatial Lyapunov exponent $\lambda$: The results of direct numerical simulations agree with the theory for all Ku if $\kappa\ll 1$ but clearly show deviations for larger values of $\kappa$. For increasing Ku, the theoretical formula (19) is approached.

In conclusion, the analytical formulas for $\lambda$ and $\mathcal{J}$, Eqs. (17) and (19), agree excellently with the results direct numerical simulations performed at finite values of Ku and St, if the inertia parameter $\kappa$ is small.

Finally we discuss the origin of the deviations between the analytical formulae [Eqs. (17) and (19)] and our numerical results when Ku is large but finite, and $\kappa$ is not small. These deviations originate from two different effects: First, the migration time between trapping regions is not negligible in general – which changes the gradient distribution so that Eq. (7) fails. Second, outside the persistent-flow limit the fluid velocity gradient changes during the relaxation of the particle velocity, so that Eq. (12) breaks down.

Our theory uses the gradient distribution $P(A_{n}=a)$ as an input into Eqs. (13) so that the importance of the first effect can be tested by modifying $P(A_{n}=a)$ appropriately. While $P(A_{n}=a)$ is given by Eq. (7) in the persistent-flow limit, it is generally determined by a complicated interplay between the dynamics of the inertial particles and the fluid motion. For that reason, so we do not have a closed-form expression for $P(A_{n}=a)$ for finite Ku and for values of $\kappa$ that are not small. Instead, we compute $P(A_{n}=a)$ from our direct numerical simulations by measuring $A_{n}$ along particle trajectories. Using the numerically computed gradient distribution as an input into Eqs. (13), we then compute the Lyapunov exponent and the rate of caustic formation with the help of the first expressions in Eqs. (16) and (18). We denote the Lyapunov exponent computed in this way by $\lambda_{A}$, and the rate of caustic formation by $\mathcal{J}_{A}$. Note that these quantities are still calculated by keeping the gradients constant on the time scales of the particle dynamics, only the gradient distribution is adjusted. Comparing $\lambda_{A}$ and $\mathcal{J}_{A}$ to $\lambda$ and $\mathcal{J}$ obtained from our direct numerical simulations thus allows us to quantify the relevance of the finite migration time between trapping regions.

Our simulation results for $P(A_{n}=a)$ are shown in Fig. 4. We observe convergence towards Eq. (7) for large Ku, but positive gradients have non-zero probability for all finite Ku. The probability for positive gradients decreases rapidly, however, as the persistent-flow limit is approached. For $\kappa=0.12$ and $\text{Ku}=10^{3}$, for example, the probability of observing positive gradients is very small, $P(A_{n}>0)\approx 4\cdot 10^{-3}$.

The results for $\lambda_{A}$ and $\mathcal{J}_{A}$ are shown as the diamonds, triangles and squares in Fig. 3(a) and (b), respectively (see Figure caption for details). We observe that $\lambda_{A}$ and $\mathcal{J}_{A}$ agree qualitatively with the results for $\lambda$ and $\mathcal{J}$ from our direct numerical simulations. Surprisingly, the path-coalescence transition for $\lambda$ is reproduced quite accurately by $\lambda_{A}$. This means that deviations from the persistent-flow limit seen in Fig. 2 at large but finite Ku and small but finite $\kappa$ are caused predominantly by the fact that the distribution $P(A_{n}=a)$ of gradients changes because the migration between traps must be taken into account. Fig. 3 shows that the persistent approximation works very well for our model, deviations occur only at large values of $\kappa$. This is in contrast to white-noise models where the Lyapunov exponent and the rate of caustic formation are strongly influenced by the fluctuations of the fluid-velocity gradients experienced in the past Gustavsson and Mehlig (2016).

IV Conclusions and discussion

We have analysed a one-dimensional, statistical model for heavy particles in turbulence in the limit of large flow persistence and weak particle inertia. Because the one-dimensional fluid-velocity field is compressible, it exhibits long-lived trapping regions where particles tend to accumulate. As the particles spend most of their time close to these traps, we have shown that the fluid-velocity gradient statistics in these regions determine the spatial Lyapunov exponent $\lambda$ and rate of caustic formation $\mathcal{J}$. We find that the rate of caustic formation is an increasing function of the inertia parameter $\kappa=\sqrt{3}\text{Ku}\,\text{St}$ in the persistent-flow limit, with exponential activation. The spatial Lyapunov exponent, by contrast, decreases monotonously as $\kappa$ increases and remains always negative. This indicates strong spatial clustering of the particles on point-like sets.

Our analytical results for $\lambda$ and $\mathcal{J}$ agree with those of direct numerical simulations of our model at large Ku as long as $\kappa$ is small. For larger values of $\kappa$ we observe deviations from the analytical results derived in the persistent-flow limit. Our analysis shows that these deviations result from the fact that migration between traps must be taken into account. This changes the positive-gradient tail of the distribution of particle-velocity gradients. We have not yet found a reliable analytical approximation for these tails.

The persistent approximation remains quite accurate even for large values of the inertia parameter $\kappa$. This means that the Lyapunov exponent and the rate of caustic formation are essentially determined by instantaneous, local properties of the flow. This is in contrast to the white-noise limit, where the history of fluid-velocity gradients experienced in the past matters Gustavsson and Mehlig (2016).

The existence of the long-lived trapping regions in the flow is an important ingredient for the one-dimensional persistent-flow model. In statistically homogeneous and isotropic flows such trapping regions exist only if the flow is compressible. Therefore, we expect particles in higher dimensional flows to exhibit a similar behaviour to that described here only if the flow is persistent, sufficiently compressible, and statistically homogeneous and isotropic.

Let us briefly discuss the relevance of the present results for incompressible, higher-dimensional flows, and for incompressible, fully developed turbulence in particular. As mentioned in Section II, the Kolmogorov time $\tau_{\text{K}}$ is analogous to the advection time scale $\tau_{\text{a}}$ in our model. However, in incompressible turbulence, and in incompressible flows in general, there are no trapping regions where advected particles stay for a long time. This means that both advected and heavy particles typically travel long distances, experiencing different fluid-velocity gradients along their way. Hence, the persistent-flow model does not directly apply to heavy particles in incompressible turbulence.

There might, however, be certain situations in turbulent flows where the combined dynamics of the fluid and the particles leads to trapping of particle trajectories. In these situations it may still be a good approximation to treat the fluid-velocity gradients at the particle position as persistent compared to the particle dynamics. Examples could be bubbles trapped in long-lived turbulent vortex-tubes by Tchen’s force Tchen (1947), or heavy particles in coherent turbulent structures with large fluid-velocity gradients, but small mean flow.

V Acknowledgements

We thank Kristian Gustavsson for discussions and for his comments on the manuscript. This work was supported by Vetenskapsrådet (Grant No. 2017-3865), Formas (Grant No. 2014-585), and by the grant ‘Bottlenecks for particle growth in turbulent aerosols’ from the Knut and Alice Wallenberg Foundation, Dnr. KAW 2014.0048. The numerical computations were obtained using resources provided by C3SE and SNIC.

Appendix A Calculation of observables from Eqs.(6a)

Instead of solving the dynamics for $Z$ as done in the main text, we show here alternative derivations for Eqs. (17) and (19) using the dynamics for $X$ instead. To this end, we first write Eqs. (6a) as

[TABLE]

Because $A_{n}$ is constant, this equation can be readily solved for $X$. We obtain

[TABLE]

where $\lambda_{\pm}=-1/2\pm\sqrt{1/4+A_{n}}$ and $C_{1}$ and $C_{2}$ are constants that depend on the initial conditions. Because $\Re[\lambda_{+}]\geq\Re[\lambda_{-}]$ for all realisations $A_{n}$, either the first term in Eq. (23) dominates, or the two terms are of the same order. We may therefore assume for simplicity that $C_{2}=0$. From the remaining exponent $\lambda_{+}$ we can now compute both the Lyapunov exponent $\lambda$ and the rate of caustic formation $\mathcal{J}$. The former is obtained by averaging the real part of $\lambda_{+}$, $\Re[\lambda_{+}]$, over the realisations of $A_{n}$. We have

[TABLE]

This expression is equal to Eq. (16), if we substitute the gradient distribution (7) for $P(A_{n}=a)$ and make the change of variables a\to z(a)=z^{*}=\lambda_{+}\big{|}_{A_{n}=a}.

The rate of caustic formation $\mathcal{J}$ is obtained by the rate at which $|X|\to 0$ while $|\dot{X}|>0$, averaged over the realisations of $A_{n}$. For fixed $A_{n}$, $X$ goes to zero at rate $\Im[\lambda_{+}]/\pi$, because $X$ passes zero twice in time intervals of length $2\pi/\Im[\lambda_{+}]$. Thus, we obtain for $\mathcal{J}$,

[TABLE]

which, again, is equal to Eq. (18) after substituting the gradient distribution (7) for $P(A_{n}=a)$, changing variables to $a\to-a-1/4$, and performing an integration by parts.