Effective rates in dilute reaction-advection systems for the annihilation process A+A\to\varnothing

# Effective rates in dilute reaction-advection systems for the annihilation process A+A→∅

G. Krstulovic Giorgio Krstulovic Jérémie Bec Laboratoire Lagrange, UMR7293, Université de Nice Sophia-Antipolis, CNRS, Observatoire de la Côte d’Azur, BP 4229, 06304 Nice Cedex 4, France Massimo Cencini CNR, Istituto dei Sistemi Complessi, Via dei Taurini 19, Roma, Italy    M. Cencini Giorgio Krstulovic Jérémie Bec Laboratoire Lagrange, UMR7293, Université de Nice Sophia-Antipolis, CNRS, Observatoire de la Côte d’Azur, BP 4229, 06304 Nice Cedex 4, France Massimo Cencini CNR, Istituto dei Sistemi Complessi, Via dei Taurini 19, Roma, Italy    J. Bec Giorgio Krstulovic Jérémie Bec Laboratoire Lagrange, UMR7293, Université de Nice Sophia-Antipolis, CNRS, Observatoire de la Côte d’Azur, BP 4229, 06304 Nice Cedex 4, France Massimo Cencini CNR, Istituto dei Sistemi Complessi, Via dei Taurini 19, Roma, Italy
July 4, 2019
###### Abstract

A dilute system of reacting particles transported by fluid flows is considered. The particles react as with a given rate when they are within a finite radius of interaction. The system is described in terms of the joint -point number spatial density that it is shown to obey a hierarchy of transport equations. An analytic solution is obtained in the dilute or, which is equivalent, the long-time limit by using a Lagrangian approach where statistical averages are performed along non-reacting trajectories. In this limit, it is shown that the moments of the number of particles have an exponential decay rather than the algebraic prediction of standard mean-field approaches. The effective reaction rate is then related to Lagrangian pair statistics by a large-deviation principle. A phenomenological model is introduced to study the qualitative behavior of the effective rate as a function of the interaction length, the degree of chaoticity of the dynamics and the compressibility of the carrier flow. Exact computations, obtained via a Feynman–Kac approach, in a smooth, compressible, random delta-correlated-in-time Gaussian velocity field support the proposed heuristic approach.

###### Keywords:
Chemical reactions; dilute media; transport; large deviations; Kraichnan ensemble
journal: J Stat Phys

## 1 Introduction

Many natural and industrial processes involve the reaction or collision of diffusing species transported by an outer flow. Such systems are typically modeled in terms of the reaction-diffusion-advection equation for the density field . In the simple case of the pair-annihilation reaction , where two molecules react together to become inert, this kinetic equation takes the form

 ∂tρ+∇⋅(ρ→v)=−Γρ2+κ∇2ρ, (1)

where is the diffusion constant, the velocity of the outer flow, and the reaction rate. An important aspect withheld in such an approach is the assumed relation between the microscopic stochastic rate, that is the probability that two given individual molecules react, and the mesoscopic reaction propensity, that is the number of reactions per unit time and volume written here as . To derive (1) one must assume two properties to be satisfied at the coarse-graining scale from which the hydrodynamic limit is taken vanKampen1992 (). First, volumes at the coarse-graining scale have to contain sufficiently many particles, in order to safely disregard finite-number fluctuations. Second, each particle within this volume must have an equal probability to react with all the others — well-mixing hypothesis. For this second condition to be satisfied, the coarse-graining scale has to be sufficiently small, i.e. smaller than the interaction distance between particle pairs.

It is clear that in situations where the reactants are very dilute, the two conditions underlying (1) might not be simultaneously satisfied. In particular, fluctuations due to a finite number of reactants might be so important to invalidate the mean-field assumptions leading to (1). Much work has been devoted to model and study such fluctuations in situations where transport is negligible and diffusion dominates. It was shown that finite-number effects can be taken into account by adding an imaginary noise in the reaction-diffusion equation. The time evolution of statistical quantities is then obtained by averaging with respect to this noise (see mattis1998uses (); Tauber2005 () for reviews). Two equivalent approaches were independently developed using either the Poisson representation gardiner1977poisson () or a field-theoretical description doi1976stochastic (); peliti1985path (). In both cases, it is assumed that a reaction takes place with a given rate when two particles are on the same node of a lattice. Diffusion is then modeled by a random walk. After that, the hydrodynamic limit is obtained by performing the limit of vanishing lattice spacing. Such approaches, which were developed for diffusion-limited reactions, cannot be straightforwardly extended to cases where transport cannot be neglected.

In many natural situations, the dynamics of very dilute reacting species is dominated by (possibly compressible) advection and their diffusion is almost negligible. This is the case, e.g., of dust grains growing by accretion to form planets in the early solar system johansen2007rapid () or of phytoplankton confined in a two-dimensional ocean layers Perlekar2010 (). Also, motile microorganisms in aquatic environments can detach from the fluid trajectories thanks to swimming and behave as tracers in a weakly-compressible flow Neufeld2007 (). Because of the large sizes of the transported species, molecular diffusion is negligible in all aforementioned examples. Another instance concerns the initiation of rain in warm clouds. Droplets, which are typically occupying a volume fraction of the order of , grow by coalescence to form raindrops. Because of their finite size and mass, such droplets have inertia and their dynamics decouples from the flow. When they are sufficiently small, droplets behave as tracers in an effective compressible flow that depends on the local turbulence of the carrier airflow Maxey2006 (). This leads to strong and violent fluctuations both in their local spatial concentration shaw2003particle () and in the rate at which they collide falkovich2002acceleration (). Classical approaches consist in using the (mean-field) Smoluchowski coagulation equation, which is a generalization of the kinetic equation (1) to the case of a full population with a size distribution pruppacher1998microphysics (). Many questions remain open on the statistical effects of turbulence on the timescales at which rain is formed devenish2012droplet (). In addition very little is known on the finite-number fluctuations induced by diluteness that might also play an important role in assessing the average growth rate of droplets. The aim of this work is to develop a general framework which can be used to address such issues. For that, we neglect diffusion and focus on particles transported by generic compressible flows. Also, we focus on the annihilation reaction because we expect that this simple model will capture the essence of finite-number effects in advection-reaction systems.

It is clear that, at long times, the particle number density decreases, so that this long-term asymptotics is very likely affected by finite-number effects. Kinetic approaches using (1) predict that, ultimately, the average number of particles decays algebraically as  chertkov2003boundary (). This behavior, which occurs at times much longer than the fluid velocity correlation time, can be explained in terms of an effective eddy diffusion. The flow compressibility might decrease the effective diffusivity vergassola1997scalar (), but the latter remains anyway positive goudon2004homogenization (). Hence, at long times, the spatial fluctuations of the concentration are smoothed out and a closed equation can be written for the spatial average of . Here, we find that finite-number effects enhance the long-time decay of the average number of particles. For that we consider discrete particles that are tracers of the compressible fluid flow and that annihilate with a rate when they are separated by less than an interaction distance . We show that the average number of particles does not decrease algebraically but rather exponentially as . This law pertains to the statistics of the relative motion between two reactants and the exponential decay rate depends on both the microscopic rate and the flow statistical properties in a non-trivial manner. The main result of this work is the introduction of a novel Lagrangian approach in terms of non-reacting particle trajectories. We exploit these ideas to express using a large-deviation principle for the time that two tracers spend at a distance below the radius of interaction .

The paper is organized as follows. Section 2 contains basic definitions and set the general framework of this work. The -point number-density field is introduced and is shown to obey a hierarchy of transport equations. When integrated over space, these fields correspond to the factorial moments of the number of particles that are present in the system. We then briefly discuss the zero-dimensional case and its relationship with standard studies of finite-number effects in well-mixed settings. In section 3, the hierarchy for the -point number-density is solved by using a Lagrangian approach that consists in following the flow characteristics. At long times the particle number moments are shown to decay exponentially with a rate that does not depend on their order. An analytic expression for is given in terms of the Lagrangian statistics of (non-reacting) tracer trajectories. In section 4 we focus on the relevant case of two particles. The exponential decay rate is then related to the time spent by the pair at a distance less than the interaction radius and is expressed in terms of the large fluctuations of the latter. Asymptotic arguments and a phenomenological picture valid for generic flows are then developed to relate the decay rate to the microscopic rate . Section 5 contains an application to very dilute suspensions in smooth compressible Gaussian delta-correlated-in-time flows (belonging to the so-called Kraichnan ensemble; see, e.g., FalkovichG.2001 ()) for which the effective reaction rate can be analytically derived. Finally, section 6 presents some concluding remarks and open questions.

## 2 Master equations in the presence of advection

### 2.1 Settings

We consider a set of particles whose positions obey the equations

 d→Xidt=→v(→Xi,t), (2)

where is a prescribed -dimensional differentiable velocity field with given isotropic, stationary, and homogeneous statistics. The velocity field is assumed to have a finite correlation time and might be compressible with compressibility

 ℘=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯(∇⋅→v)2 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯tr% (∇→vT∇→v) , (3)

where the over line designates the average with respect to the velocity field realizations, and . The extremal values and correspond to incompressible and purely potential velocity fields, respectively. Whenever , the dynamical system defined by (2) is dissipative. Also, to maintain sufficiently generic settings, the dynamics is assumed ergodic and chaotic. Typically, when the flow is defined on a compact set, the trajectories generated by (2) will concentrate on a dynamically evolving strange attractor, while in the incompressible limit they remain uniformly distributed (see, e.g., FalkovichG.2001 () for more details on particle transport in compressible and incompressible flows). For the sake of simplicity, in this paper we will focus on the dynamics of tracers given by (2). However, the results can be straightforwardly extended to a general dynamics given by the Newton equation as that ruling, for instance, the evolution of inertial particles Maxey2006 ().

The dynamics (2) is supplemented by binary reactions between the particles. When the particles labeled and might react and annihilate (or become inert). This happens with a rate . As a consequence, the number of particles in the domain will decrease from its initial value in the course of time. We suppose that the radius of interaction is smaller than the scale at which varies. For example, in turbulent flows, this corresponds to assuming that is smaller than the Kolmogorov length scale . Moreover, in writing (2) we have neglected particle diffusion. This is justified whenever the interaction distance is larger than the Batchelor length scale , where designates the Schmidt number — that is the ratio between the fluid kinematic viscosity and the particle diffusivity . For applications, our settings thus reduce to very large values of the Schmidt number. We notice that the Schmidt number can be of the order of thousands or higher in organic mixtures, biological fluids and generically for particulate matter. Estimates based on Stokes–Einstein relation lead for a micro-meter sized spherical particle in air and in water.

Before discussing the master equations which describe the system, it is worth underlining that our settings involve different sources of stochasticity. First, the dynamics (2) has to be supplemented by initial conditions on particle positions that we choose randomly. A second source of randomness comes from the intrinsic stochasticity of the reaction process. Finally, the fluid flow is stochastic with prescribed statistics. Consequently, ensemble averages corresponding to these different sources of randomness must be considered; the corresponding notations will be introduced when needed.

### 2.2 Master equations for n-point number densities

To fully describe all interactions between the particles in the considered system, all relative positions are simultaneously needed. It is then natural to define for ( being the initial number of particles) the joint -point number density

 Fn(→x1,…,→xn,t)=⟨∑i1≠…≠inδ(→Xi1(t)−→x1)⋯δ(→Xin(t)−→xn)⟩μ, (4)

where the brackets stands for the ensemble average with respect to reactions only and, by convention, the sum in (4) is zero when . The particle indices vary between and the number of particles in each realization of the reactions. The sum is not ordered so that, by construction, the spatial averages of the ’s are the factorial moments of the total number of particles, namely

 F1(t) = (5) F2(t) = ∬F2(→x1,→x2,t)ddx1ddx2=⟨N(t)(N(t)−1)⟩μ (6) Fn(t) = ∫⋯∫Fn(→x1,…,→xn,t)ddx1⋯ddxn=⟨N(t)!(N(t)−n)!⟩μ. (7)

Hence, is the average number of -uplets among the particles still present at time .

Taking into account the transport by the flow and the reactions, the joint -point densities (4) satisfy the following hierarchy of equations

 ∂Fn∂t+n∑i=1∇→xi⋅[→v(→xi,t)Fn] = −μFn∑i

where is the -dimensional divergence with respect to the spatial variable and the Heaviside function. The transport term comes from the identity

 ddtδ(→Xik(t)−→xk)=−∇→xk⋅[δ(→Xik(t)−→xk)→v(→xk,t)]. (9)

The two terms in the right-hand side are sinks accounting for the reactions: the first includes all binary reactions between any close-enough pair among the considered particles; the second counts the reactions of any of these particles with another one at a distance smaller than .

The hierarchy (8) for the -point joint number density is exact and contains all information about the dynamics and statistics of the system. Moreover, it is closed at the order and can be integrated by the method of characteristics (see Sect. 3). In the next subsection we relate this description in terms of factorial moments to standard descriptions used in the statistical physics of well-mixed chemical reactions.

### 2.3 The zero-dimensional case

When the hypothesis of well-mixedness is satisfied, as when for instance the interaction radius is larger than the system size, transport does not play any role in the temporal evolution of the factorial moments (4). In this zero-dimensional limit, the hierarchy of transport equations (8) becomes a simpler hierarchy of ordinary differential equations, i.e.

 (10)

which can be also derived from the standard master equation describing the reaction . Denoting with the probability of having particles at time , we have vanKampen1992 ()

 dPNdt=μ2(N+2)(N+1)PN+2−μ2N(N−1)PN. (11)

By definition, the factorial moments are expressed in terms of as

 Fn(t)=∞∑N=0N!(N−n)!PN(t). (12)

Taking the time derivative of and using (11) yields to (10).

Equation (11) is generally the starting point for describing the statistics of chemical reactions. In the limit of large number of particles, the mean-field approximation is often used to write from(11) a closed equation for the mean number of particles. This mean-field equation, obtained by replacing by in (10) for , predicts a algebraic decay of the number of particles. Clearly, the mean-field approximation will breakdown at long times when the number of particles is smaller. Many works have been devoted to describe fluctuations and deviations from mean field due to a finite number of particles CardyImaginaryNoise (); Deloubriere2002 (); Tauber2005 (). In particular, it is possible to show that solving equation (10) is equivalent to computing the moments of the solution of a stochastic differential equation with a pure imaginary noise CardyImaginaryNoise (); Gardiner_1985 (). However, this approach cannot be directly applied when the interaction radius is finite and transport is present. Reactions and spatial location of particles cannot be treated separately as it is usually done in reaction-diffusion systems, where the well-mixed hypothesis is assumed at the scale of the lattice spacing. The hierarchy of transport equations (8) naturally couples the spatial distribution of particles and the reactions. The aim of the next section is to provide a formal solution of (8) in terms of Lagrangian trajectories (or tracers) of the carrier flow. These trajectories are completely defined by carrier flow and do not depend on the reactions.

## 3 A Lagrangian approach to particle reaction kinetics

### 3.1 Finite-number closure and recursive solutions

We now consider a finite interaction distance so that the presence of transport by the flow must be explicitly accounted for. Without loss of generality, we can assume that at the total number of particles is fixed and equal to . As for all , the hierarchy (8) is closed at the order and can be integrated.

 ∂FN0∂t+N0∑i=1∇→xi⋅[FN0→v(→xi,t)]=−μFN0∑i

which can be formally solved in term of Lagrangian trajectories — that is using the method of characteristics. Denoting by the solution to the differential equation

 ddt→Y(t;→y)=→v(→Y(t;→y),t),→Y(0;→y)=→y, (14)

it is straightforward to check that the solution of (13) is given by

 FN0(→x1,…,→xN0,t)= =∫F0N0e−μ∑i

with . We have introduced here the so-called Lagrangian average , which is an average over all -uplets of non-reacting trajectories such that their location at time is . Such an average is often used in compressible transport gawedzki2000phase (). We can now insert the expression for in the equation for obtaining

 ∂FN0−1∂t+N0−1∑i=1∇→xi⋅[FN0−1→v(→xi,t)]=−μFN0−1N0−1∑i

where

 Λn(t)=μ∑i

 (18) −μ⟨F0N0e−ΛN0−1(t)N0−1∑i=1∫t0θ(a−|→Y(s;→yi)−→Y(s;→yN0)|)eΛN0−1(s)−ΛN0(s)ds⟩N0.

The general solution for can be obtained reintroducing (18) in (8) solving for and iterating until the order .

In writing (14) we have implicitly assumed the existence and uniqueness of the Lagrangian trajectories and each trajectory is labeled by its initial position. This assumption is guaranteed as we have considered a smooth velocity field. For a fully turbulent flow and if is much larger than the Kolmogorov length scale, the velocity field is no longer smooth but only Hölder continuous with an exponent smaller than one. In such a case the solutions to (14) are not unique and the deterministic Lagrangian flow breaks down. However, a statistical description of Lagrangian trajectories is still possible gawedzki2000phase ().

We conclude this subsection by noticing that this new approach in terms of Lagrangian trajectories allows us to directly apply known results from passive (non-reacting) transport (see FalkovichG.2001 () for a review) to address the problem of advection-reaction. These theoretical considerations will be mostly used in Sect. 4.

### 3.2 Long times and very dilute closures

When considering very dilute systems and small interaction radii, the hierarchy of transport equations can be closed at each order as the second term of the right hand side of Eq. (8) can be neglected with respect to the first one. In particular, the closure becomes exact in the long time limit when the mean number of particles depends only on the statistics of the relative distance of one pair of Lagrangian trajectories.

The long time behavior of the joint -point density can be directly obtained from equation (18). Noticing that , at the leading order, the joint -point density is given for by

 Fn(→x1,…,→xn,t)=⟨e−Λn(t)⟩n. (19)

Integrating over all spatial variables, for , we obtain

 Fn(t)=⟨F0ne−μn∑i

where now the average is over all possible -uplets of tracers. Finally, averaging the factorial moments over the realizations of the velocity field, we can define the effective reaction rate as

 γn=−limt→∞1tln¯¯¯¯¯¯¯¯¯¯¯¯Fn(t)¯¯¯¯¯¯¯¯¯¯¯¯¯Fn(0). (21)

Note that the average over different realizations of the carrier flow, denoted with an overline, is performed before taking the logarithm. This important point will be further discussed in Sect. 4.1.

In the zero-dimensional limit, by definition the Heaviside function inside equation (20) is equal to one, and thus the effective rate is simply given by , as confirmed by numerical data (not shown). When the system is not well-mixed at all scales, the effective rate (21) cannot be computed for generic flows as multi-time and space correlations of Lagrangian trajectories are involved. However, an important role is played by the case that only depends on the two-point motion. Indeed, it is easy to check from equation (20) that is an increasing function of . In addition, the moments of the number of particles can be expressed as a linear combination of the factorial moments. The rate is thus always the leading exponent and for all we have

 ⟨Np⟩∼e−γ2t. (22)

This exponential decrease is expected to be valid at long times regardless of the initial number of particles considered. This behavior is evidenced in Fig. 1 where the temporal evolution of the moments is displayed for , and different orders . For , the long-time departure from the mean-field prediction is clear from the inset.

Independently of the order of the moment considered, the long-time dynamics depends only on the statistics of the relative distance of one pair of tracers, which motivates the next section that is devoted to the study of the two-point motion.

## 4 Two-point motion

### 4.1 Effective rate

As the long-time behavior of the system is dominated by the two-particle dynamics, it is natural to close the hierarchy (8) at the second order and to look at the mean numbers of particles and of pairs:

 n1(t)≡¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯⟨N(t)⟩μ=¯¯¯¯¯¯¯¯¯¯¯F1(t) , n2(t)≡12¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯⟨N(t)[N(t)−1]⟩μ=12¯¯¯¯¯¯¯¯¯¯¯F2(t). (23)

The overbars designate the average with respect to the realizations of the velocity field , which is assumed statistically homogeneous in space so that the overbars include the average over the particle initial positions. Both quantities, and , are expected to exponentially decrease with the same rate , which here and in the following will be denoted by . An effective equation describing such a behavior is simply given by

 dn1dt=−γn2 , dn2dt=−γn2. (24)

The effective rate depends on the microscopic rate , the interaction radius and the statistics of the carrier flow through equations (20) and (21). However, these equations do not allow us to readily obtain an explicit expression for . The dependence of such an expression on the parameters of the system can be obtained by using standard tools of dynamical system and statistical physics.

The effective rate , defined by equation (21), depends on the particles trajectories only through and thus is determined by the two-particle statistics. Consequently, without loss of generality, we will focus on a system initially having particles. For such a system, and the average number of pairs reads

 n2(t)=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯e−μ∫t0θ(a−R(s))ds, (25)

where and the average is over all realizations of the carrier velocity field . The argument of the exponential in Eq. (25) contains a time integral that, for ergodic two-point dynamics, self-averages, i.e.

 limt→∞1t∫t0θ(a−R(s))ds=P<2(a)≡P{R(s)

where is a shorthand notation for probability. This observation leads to the following naive prediction

 γnaive=μP<2(a). (27)

This estimate is consistent with what could be expected on a heuristic basis, however it is valid only when considering a fixed realization of the carrier flow. This kind of setting corresponds to the so-called quenched disorder in statistical mechanics. Indeed, let us fix a realization of the carrier flow and the initial positions and of the pair of reacting particles. If we now denote by the variable which is equal to when the two particles have not reacted and otherwise, we clearly have that

 limt→∞1tln⟨N(t)[N(t)−1]⟩μ=−limt→∞μt∫t0θ(a−R(s))ds=−μP<2(a).

Notice that, here, the average is performed solely with respect to reactions. The quenched effective rate is thus given by (27). The average over fluid realizations can be subsequently performed but it does not modify the result (27) as is a self-averaging quantity.

The situation is quite different when considering an annealed disorder, that is when the fluctuations of the carrier flow are taken into account to compute the number of pairs, according to the definition (25). In this case, as is not a self-averaging quantity, fluctuations of the carrier flow affect the long-time behavior and, hence, is not simply given by (27). In most physical situations, fluctuations of the carrier flow and randomness of the initial particle positions cannot be disregarded, so that annealed averages must be performed to define the effective rate (21). As discussed in the next subsections, in this case, the effective rate can be related to the large deviations of the fraction of time spent by two particles within a distance less than .

### 4.2 Large deviations of the reaction time

Determining the long-time asymptotics of requires to understand the behavior of the finite-time average

 Θ=1t∫t0θ(a−R(s))ds. (28)

Under ergodicity hypotheses on the two-point dynamics, this self-averaging quantity converges at long times to . When the dynamics is sufficiently mixing, the deviations of from its average are described by a large-deviation principle Ellis85 (). The probability density function (PDF) thus has the asymptotic behavior

 −limt→∞1tlnp(Θ)=H(Θ), (29)

where is a convex positive rate function, which attains its minimum, equal to zero, at . Plugging this form in (25) we obtain that

 n2(t)∝∫∞0e−t(μΘ+H(Θ))dΘ. (30)

A saddle-point estimate of the integral shows that the effective rate is given by the Legendre transform of the rate function :

 γ=infΘ≥0[μΘ+H(Θ)]. (31)

The convexity of implies that is a non-decreasing function of . For illustrative purposes, in figure 2 (A) we show the rate function for the Kraichnan flow (a delta-correlated in time velocity field, see Sec. 5 for numerical settings and throughout discussion of this kind of flows).

When the reaction rate is small compared to the scale of variation of , the infimum in (31) is attained very close to the minimum of at . There, as a result of the central limit theorem (CLT), the rate function is approximately quadratic (represented by the black solid line in figure 2 (A)) and the infimum in (31) is attained for , so that

 γ≃μP<2(a)−μ2σ22, (32)

where is the variance of the fraction of time spent by the pair at a distance less than . This behavior is apparent figure 2 (B) where the effective rate is displayed for the same flow. The expression (32) is valid as long as is close enough to to neglect the sub-leading terms in the quadratic approximation for , i.e. in the domain of validity of the CLT. Notice that (32) requires . Observe also that the naive estimate (27) for a quenched velocity field only captures the linear behavior of close to .

The asymptotic behavior of at large values of relates to that of the rate function at small values of the fraction of time spent below , indeed, from (31) we have

 limμ→∞γ=limΘ→0H(Θ). (33)

The events leading to small ’s correspond to those in which the particle-pair separation remains larger than for a very long time. More precisely, the value of at relates to the distribution of the random time needed by pairs that are initially far apart (e.g. at separations of the order of the system size) to approach each other at a distance less than . We can indeed write

 P(Θ<ε)=P(Θ<ε|Thit>t)P(Thit>t)+P(Θ<ε|Thit

When , the two particles have never been at a distance less than , which implies for any . Conversely, when , the fraction of time is finite and for . Hence, it follows that

 H(0)=−limt→∞1tlnP(Thit>t). (35)

Different behaviors of the effective rate are thus expected depending on the tail of the distribution of . For the proposed large-deviation approach to be valid, one expects to decrease at least as fast as an exponential. Otherwise, one would obtain , violating the convexity of . When decreases faster than an exponential, , and the effective rate for . The way it diverges cannot be obtained from the distribution of only, as it depends on the functional form of close to . In the intermediate case, when decreases as an exponential, is finite and the effective rate approaches a finite value when . Depending on whether is finite or not, saturates to for or approaches it asymptotically. Note that an exponential behavior for the tail of the distribution of is not a particular case but is actually expected to be a quite generic circumstance. Indeed, when considering times much longer than all relevant timescales of the flow, memory loss implies that the hittings of form a Poisson process, and the hitting time is an exponential random variable. This is directly corroborated when considering the Kraichnan flow (data not shown). For an exponential distribution of the hitting time we simply obtain

 γ∞=1¯¯¯¯¯¯¯¯Thit. (36)

This asymptotic behavior is also observed in figure 2 (B) for large values of .

In the next subsection we introduce and study a phenomenological model for which the rate function can be explicitly computed.

### 4.3 A phenomenological model for reactions in bounded chaotic flows

The essence of chaotic flows is the competition between stretching and folding. Stretching results from the fact that close trajectories separate exponentially at a rate given by the largest Lyapunov exponent . Folding, which is generally due to boundary conditions, prevents the inter-particle distance to grow indefinitely and is necessary for the system to converge to a statistical steady state. Qualitatively, the time evolution of the distance between two particles resembles the trajectory shown in Fig. 3. The particle pair spends long times at a distance of the order of the size of the domain and performs rare excursions to very small separations.

When interested in particles reacting within a distance , it is quite natural to decompose the inter-particle dynamics in a sequence of consecutive phases, each corresponding to a given approaching event and separated by returns to . Starting from , a first time subinterval is required to hit for the first time. Then the particle separation grows to touch again . The time interval is thus decomposed in random subintervals of length , where is the hitting time from to and the time needed to reach for the first time starting from (see figure 3). Because of the memory loss occurring when the pairs are at a distance , these sub-intervals can be considered independent and identically distributed. The first substage corresponds to many recurrences to during time lengths that are independent from each other because of the imposed boundary condition. The hitting time is dominated by the sum of these recurrence times and this leads to assume that has an exponential distribution of mean . Concerning the escape time , it is clear that when , its behavior is dominated by the long-time convergence to the Lyapunov exponent. We can thus assume that it is deterministic, namely . In addition, in this limit, we always have and thus

 t≈nrea∑n=1[Tnhit+Tnesc]≈nrea∑n=1Tnhit. (37)

According to the discussion of the previous subsection, the times are exponentially distributed and is a Poisson random variable of mean .

During the -th interval of length , the particles spend a time below a distance . The total fraction of time spent at can thus be written as

 Θ=1tnrea∑n=1αnTnesc≈nreat¯¯¯¯αTesc, (38)

where designates the average fraction of time spent below during . We have assumed in this model that the fluctuations of do not play an important role. Notice that under these assumptions, . We see from equation (38) that is proportional to . Its probability density is thus given by the Poisson distribution. Using the Stirling formula when , we obtain

 H(Θ)=1¯¯¯¯¯¯¯¯Thit−Θ¯¯¯¯αTesc(1+ln[¯¯¯¯αTescΘ¯¯¯¯¯¯¯¯Thit]). (39)

It is easy to check that given by (39) is a convex function that vanishes at its minimum reached at , as expected. Also, one sees that and , where the prime denotes the derivative with respect to . In the light of the considerations of the previous section, this implies that the effective reaction rate asymptotically converges to when . With the explicit form (39) for the rate function we actually obtain

 γ=1¯¯¯¯¯¯¯¯Thit[1−e−P<2(a)¯¯¯¯¯¯¯¯Thitμ]. (40)

Note that this expression also reproduces the small- asymptotics (32) with .

The explicit form (40) relies on the strong simplifications of the particle dynamics that we have made (cf. Fig. 3). However, this model is particularly useful to understand how depends on the reaction radius , the flow compressibility, and its chaoticity that will affect both and . The time can be seen as the recurrence time from to . Therefore, by the Kac Recurrence Lemma KacLemma (), it is expected that . However, the proportionality factor depends on the Lyapunov exponent and on the compressibility, and thus on the detailed properties of the carrier flow.

The model presented in this section, although phenomenological, is based on quite general assumptions, such as the existence of a finite correlation time and of at least one positive Lyapunov exponent. The precise shape of the rate function (39) can vary from one system to another when more realistic dynamics is considered. It is however reasonable to expect that the general qualitative properties of the effective rate remain essentially the same.

## 5 Applications to very dilute suspensions in compressible Kraichnan flows

To illustrate and further validate the approach developed in the previous section, we focus here on two-particle motion in a special class of random velocity fields, that is the smooth (compressible) Kraichnan ensemble gawedzki2000phase (); FalkovichG.2001 (). The main advantage of such a class of stochastic flows is that, thanks to time uncorrelation, we have analytical control on most of the quantities we need. Further it is very efficient to implement numerically so to gain high quality statistics.

### 5.1 The Kraichnan velocity ensemble

We focus on dimensions . The velocity is a homogeneous, isotropic, Gaussian spatially smooth field, with zero mean and two-point correlation function

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯vi(→0,t)vj(→r,t′)=2Dij(→r)δ(t−t′),Dij(→r)=D0δij−dij(→r)/2. (41)

For small separations the spatial part of the correlation is prescribed according to gawedzki2000phase ()

 dij(→r)=D1[(d+1−2℘)δijr2+2(℘d−1)rirj]. (42)

The flow parameters are as follows: the diffusivity controls the single particle dispersion properties; has the dimension of an inverse time and sets the intensity of velocity gradients; finally, controls the compressibility degree as defined in Sect 2.1, its extreme values correspond to incompressible and potential velocity fields, respectively. When advecting particles with the velocity field defined by Eqs. (41)-(42), we will always consider the presence of boundary conditions, which are necessary in order to ensure statistical steady state properties (see also below the description of the numerical implementation).

As for the numerical implementation we consider two possible schemes. The first, which is useful when considering many particles (this is that used for Fig. 1), consists in generating a random velocity field, with periodic boundary conditions, as the sum of independent Fourier modes with and amplitudes such that Eqs.(41)-(42) are satisfied at small scales. When interested in two-particle motion, we use a second scheme, which is much more efficient. It evolves directly the stochastic equation for the separation between the particles. We use a discretization by means of a standard Euler-Ito scheme that reads

 Ri(t+dt)=Ri(t)+√2dtVi(→R,t) for i=x,y with Vx=√D1(1+2℘)Rxηx−√D1(3−2℘)Ryηy with Vy=√D1(1+2℘)Ryηx+√D1(3−2℘)Rxηy,

which only requires to generate two independent zero-mean Gaussian random variables, . The above equations indeed prescribe that , with given by Eq.(42) for . Finally, to ensure a statistically steady dynamics, reflective boundary conditions are imposed at . In principle, this scheme can be generalized also to more than two particles FMNV99 ().

Particles evolving in such flows with the dynamics Eq. (2) define a stochastic dynamical systems of which we know explicitly the Lyapunov exponents, (), and the rate function of the stretching rates, which is quadratic (see the review FalkovichG.2001 () for a summary of the known properties of Kraichnan flows). In particular, the particle separation follows a diffusion process whose generator, restricted to homogeneous and isotropic sectors, is simply given by

 M2=D1(d−1)(1+2℘)rD2−1∂∂r(rD2+1∂∂r), (43)

where

 D2=d−4℘1+2℘=d−℘2(d+2)1+2℘. (44)

To ensure the convergence to a statistically stationary regime, the particle dynamics and the generator must be defined on a compact set; this is done by supplying boundary conditions at the domain border. The resulting stationary probability of finding two particles at a distance less than is

 P<2(a)∼(a/L)D2for a≪L (45)

meaning that defined in (44) is the correlation dimension.

It is worth noticing that for , which simply means volume preserving dynamics in the incompressible limit, and that for , which is the value of the compressibility above which the first Lyapunov exponent becomes negative. The latter property indicates that above a critical compressibility there is a transition to a regime of particles trapping with particles asymptotically collapsing onto a single point see, e.g., Ref. gawedzki2000phase (). Notice that this regime cannot be realized for .

### 5.2 Feynman–Kac approach

We now turn to a direct computation of the effective rate by using the generator (43) of the diffusion process . Note that by performing the change of variable , reduces to

 ~M2=λ1∂∂x+Δ2∂2∂2x, (46)

where is the largest Lyapunov exponent and is the variance of its finite-time fluctuations. Hence, for the Kraichnan model is a Brownian motion with drift and diffusion coefficient gawkedzki2004sticky (). Note also that the correlation dimension is expressed as

 D2=2λ1/Δ. (47)

To obtain the effective rate we need to compute the large-time dynamics of the number of pair given by equation (25). Conditioning on the initial separation between the particles, we define

 ψ(x,t)=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯[e−μ∫t0θ(log(a/L)−xs)ds∣∣xs=0=x], (48)

where is the Brownian motion with drift defined by the generator . By definition is restricted to the half-plane and we assume reflecting boundary conditions at (corresponding to the domain boundary ). The effective rate is thus given by the decay rate of . The Feynman-Kac formula donsker1975asymptotic (); oksendal1998stochastic () states that is the solution of the partial differential equation

 (49)

The boundary condition at comes from the no-flux condition expressed for the backward operator of the reflected Brownian motion with drift defined by (46) (see e.g. Gardiner_1985 ()). The effective rate is thus given by the modulus of the largest eigenvalue of the operator . The eigenvalue problem can be easily solved finding the solutions for and and then imposing that and its derivative are continuous at . Finally, imposing the boundary condition at gives the following transcendental equation for the eigenvalues :

 −(λ21−λ1√λ21+2Δ(μ+s)+2Δs)√λ21+2Δssinh⎛⎜ ⎜⎝log(a/L)√λ21+2ΔsΔ⎞⎟ ⎟⎠+ (50) (√λ21+2Δ(μ+s)−λ1)cosh⎛⎜ ⎜⎝log(a/L)√λ21+2ΔsΔ⎞⎟ ⎟⎠=0.

All solutions of equation (50) are negative and the absolute value of the largest one corresponds to the effective rate . Taking the limit in equation (50) leads to a transcendental equation for the mean value of the hitting time . In addition, solving (50) for small and using the asymptotic value of (32) directly gives and the variance of the fraction of time spent by a pair of particles at a distance less than :

 P<2(a)=(aL)D2 (51) σ2=4P<2(a)D2λ1[1−P<2(a)(1−log[P<2(a)])]. (52)

### 5.3 Effective rate for the Kraichnan flow

The Feynman-Kac calculations presented in the previous section are independent of the space dimension . This is a direct consequence of the fact that the diffusive behavior of the logarithm of the relative distance does not present any dependence on other than through and (see equation (47)). In the following, we will thus focus on the two-dimensional case, which is the simplest non-trivial example. In particular, we shall consider compressibility values to stay away from the trapping regime, which is trivial from the point of view of reaction dynamics.

As discussed at length in the previous sections, the long-time statistics of very dilute systems is dominated by the two-particle effective reaction rate . This is confirmed by numerical simulations of the compressible Kraichnan flow: as shown in Fig. 1 the decay of the moments of the number of particles is indeed asymptotically dominated by . Therefore, here, we focus on the two-point dynamics.

The effective rate as a function of for different values of the compressibility is displayed in figure 4. Both, Feynman-Kac calculations and numerics are in good agreement. Tiny discrepancies are due to the lack of numerical precision when estimating the exponential in (25) for large .

It is apparent in figure 4 that compressibility enhances reactions. The same behavior is also observed in .

The effective rate clearly exhibits two different regimes for small and large values of . For small , the observed linear behavior is given by the naive approximation (see equation (27)). This approximation remains valid as long as is much smaller than the fluctuations of the fraction of time . From the central limit theorem prediction (32), we infer that this crossover takes place for . Using (52) we thus obtain that the linear behavior on the microscopic rate, only holds for , as confirmed in figure 4. On the other hand, for large , the rates saturates to the value