Memory effects in chaotic advection of inertial particles

A systematic investigation of the effect of the history force on particle advection is carried out for both heavy and light particles. General relations are given to identify parameter regions where the history force is expected to be comparable with the Stokes drag. As an illustrative example, a paradigmatic two-dimensional flow, the von K\'arm\'an flow is taken. For small (but not extremely small) particles all investigated dynamical properties turn out to heavily depend on the presence of memory when compared to the memoryless case: the history force generates a rather non-trivial dynamics that appears to weaken (but not to suppress) inertial effects, it enhances the overall contribution of viscosity. We explore the parameter space spanned by the particle size and the density ratio, and find a weaker tendency for accumulation in attractors and for caustics formation. The Lyapunov exponent of transients becomes larger with memory. Periodic attractors are found to have a very slow, $t^{-1/2}$ type convergence towards the asymptotic form. We find that the concept of snapshot attractors is useful to understand this slow convergence: an ensemble of particles converges exponentially fast towards a snapshot attractor, which undergoes a slow shift for long times.


Abstract.
A systematic investigation of the effect of the history force on particle advection is carried out for both heavy and light particles. General relations are given to identify parameter regions where the history force is expected to be comparable with the Stokes drag. As an illustrative example, a paradigmatic two-dimensional flow, the von Kármán flow is taken. For small (but not extremely small) particles all investigated dynamical properties turn out to heavily depend on the presence of memory when compared to the memoryless case: the history force generates a rather nontrivial dynamics that appears to weaken (but not to suppress) inertial effects, it enhances the overall contribution of viscosity. We explore the parameter space spanned by the particle size and the density ratio, and find a weaker tendency for accumulation in attractors and for caustics formation. The Lyapunov exponent of transients becomes larger with memory. Periodic attractors are found to have a very slow, t −1/2 type convergence towards the asymptotic form. We find that the concept of snapshot attractors is useful to understand this slow convergence: an ensemble of particles converges exponentially fast towards a snapshot attractor, which undergoes a slow shift for long times.

Introduction
The advection of small, rigid, finite-size particles has been studied recently in several setups because it plays an important role in many environment-related phenomena ranging from atmospheric science and oceanography to astrophysics (for reviews see [1,2,3]). For example, inertial particles are important in cloud microphysics [4,5], planet formation [6], and aggregation and fragmentation processes in fluid flows [7,8,9]. A particularly interesting aspect is the advection of such particles in chaotic flows where experimental results [10,11] are available by now. Applications can be for example the forecasting of pollutant transport for homeland defense [12,7] and the localization of a toxin or biological pathogen spill (e.g. anthrax) from outbreaks in a street canyon [13].
Due to the drag acting on the particles, the Hamiltonian passive advection problem is converted into a dissipative problem which can have attractors. Correspondingly, inertial particles can have the tendency to accumulate in certain regions of the flow [14] (a phenomenon termed preferential concentration). A characteristic feature in comparison with the dynamics of ideal tracers is the appearance of caustics [4,15], i.e., the intersection of different branches of the chaotic sets in the configuration space. This is a consequence of the fact that what we see is a projection of the full pattern of the high-dimensional phase space to the configuration space of the fluid. The main relevance of the existence of caustics is that the probability of collisions of particles becomes enhanced in such regions of the flow.
The basic equation of motion for small spherical particles in a viscous fluid describable in the Stokes regime is given by the Maxey-Riley equations [16,17] ‡. Its precise form contains an integral, also called the history force, which describes the diffusion of vorticity around the particle during its full time history. The kernel of this integral is therefore proportional to the power -1/2 of the elapsed time.
The origin of such a kernel can perhaps best be seen in the problem of an infinitely large plate moving in a rather viscous fluid. If the plate starts at time zero with a sudden change in the velocity which is then kept constant, in Stokes approximation, one solves the diffusion equation for the velocity field. The solution contains a prefactor 1/ √ t. For the case of a general motion of the plate, one imagines to break up the problem in a set of jumps in the plate's velocity at different times τ . Each such jump contributes a factor to the flow velocity at time t proportional to 1/ √ t − τ multiplied by the velocity jump. In a continuous time picture, the result for the shear stress on the plate is proportional to an integral from zero to t of the plate's acceleration at time τ multiplied by the kernel 1/ √ t − τ (see , paragraph 24). For the motion of a spherical particle in a viscous flow it were Boussinesq and Basset (at the end of the 19th century) who pointed out the appearance of such a history force (for historical details see [20]). This force is therefore sometimes called the Boussinesq-Basset or the Basset ‡ It is interesting to note that Gatignol derived the same equation as Maxey and Riley in the same year. Currently, it is this equation augmented with the corrections introduced by Auton and coworkers [18] (see Eq. (1)) that is commonly called the "Maxey-Riley equation" force (we use the term history force). The equation of motion of a moving sphere in a fluid at rest was already given in , problem 7 to paragraph 24). The precise form of the equation in a moving fluid [16,17,18], i.e., Eq. (1) or (6), was, however, established only a century after the first works of Boussinesq and Basset.
The history force renders the advection equation of inertial particles to be an integro-differential equation. Because of this difficulty, the memory represented by this integral term was neglected in the majority of papers on chaotic advection of inertial particles.
In nonchaotic flows, the history force was shown to have important effects [21,22,23,24,25]. A particularly interesting example is a recent experiment which points out the determining importance of memory in the motion of microbubbles propelled by ultrasound [26,27]. In this setting the trajectories, and a surprising destabilization of these (in spite of rather strong dissipation), can be understood only when the history force is taken into account. Within a perturbative treatment restricted to the limit of very weak inertia, it has been found that the history force can suppress chaotic advection [28]. In turbulent flows, the effects of the history force have been studied in [29,30,31,32], showing that this force can have significant influence. For neutrally buoyant particles it has been suggested in [11] that there are dynamical effects beyond the history force, Faxén corrections and the lift force.
Our aim is to carry out a comprehensive investigation of long term memory effects due to strong inertia in smooth flows with chaotic advection, made possible by a recent progress in the numerical treatment of the problem [33]. Concerning the physical aspects, we are extending here the results of a short paper [34] before which, to our knowledge, no such investigations have been carried out.
The paper is organized as follows. First, we review the basic equation (Section 2), and provide estimates for parameters where memory effects are expected to be relevant (Section 3). Next we sketch the applied numerical scheme (Section 4). In Section 5, a paradigmatic two-dimensional chaotic model flow is introduced, that of the von Kármán vortex street, which shall be used as an illustrative example throughout the paper. The choice of the particle parameters is discussed in Section 6. In chaotic flows, the basic question is not so much about the deviation between trajectories with and without memory (treated briefly in Section 7), rather the deviation in statistical properties. We shall point out in Section 8 that properties like the escape rate or the lifetime statistics of particle ensembles basically differ due to the history force both for bubbles and aerosols (particles lighter and heavier than the fluid). The residence time distribution and the value of the power-law exponent of nonhyperbolic decay is also different in the two cases, as shown in Section 9. Other types of statistics, like e.g. that of the different forces acting on particles while being in the wake is investigated in Section 10. One important effect of memory is that the tendency for accumulation is much weaker than with Stokes drag only and, as a consequence, attractors occur only at rather special parameters in this open flow. When attractors nevertheless exist, the temporal convergence of individual trajectories towards them is converted from an exponential approach to a power-law approach (Section 11). It is worth, however, replacing the traditional individual trajectory picture by a global one, and monitoring the evolution of a particle ensemble. In this picture an exponential convergence is found to a so-called snapshot attractor, which then slowly moves towards the traditional attractor (Section 12). After briefly discussing also cases with chaotic attractors (Section 13), we show in Section 14 that the average Lyapunov exponents of transient chaos with memory are larger than without. The details of how to define and determine Lyapunov exponents in an integrodifferential equation, like the Maxey-Riley equation, are relegated to Appendix A. The paper is concluded with a summary of several dynamical features illustrating our main finding: memory effects change the inertial dynamics so that it moves towards that of the passive dynamics, but in a highly nontrivial way. We also return to the discussion of the relevance of the history force and argue that the size parameter is most naturally defined in terms of the shortest characteristic time scale. A refined estimate is found, valid for any flow, enabling us to also briefly discuss the relevance of memory effects in turbulent flows.

Basic equations
The equation of motion of a rigid spherical particle of radius a and mass m p in a fluid of kinematic viscosity ν, density ρ f and velocity field u(r, t) reads (including gravity) as [16,17,18] Here v ≡ dr/dt is the particle velocity, m f is the mass of the fluid excluded by the particle, g the gravitational acceleration, and and denote the full derivative along the trajectory of the particle and of the corresponding fluid element, respectively. The terms on the right-hand side of (1) are: the force exerted by the fluid on a fluid element at the location of the particle, the added mass term § describing the impulsive pressure response of the fluid, the Stokes drag, the buoyancyreduced weight, and the history force. The history force, an integral term, accounts for the viscous diffusion of vorticity from the surface of the particle along its trajectory [16]. § We note that in the original derivation by Maxey and Riley, and by Gatignol the added mass term contained the derivative du dt . Later it was found that Du Dt is a better choice as it is correct in a wider range of conditions (see [18] and references in [35]).
Equation (1) is valid if the initial particle velocity v(t 0 ) at time t 0 matches the fluid velocity. Otherwise, an initial-velocity-dependent term [36,35] should be added to the right-hand side of (1). By means of the identity (that can be derived by partial integration) we arrive at the most general form of the history force (the right-hand side of (5)) valid for any initial condition. Note that this is formally nothing but d the fractional derivative of order 1/2 of the slip velocity v − u (more precisely it is a fractional derivative of the Riemann-Liouville type [37]). A new feature of the right-hand side of (5) is that the nominator of the integrand does not longer contain a derivative. It is just the velocity difference what appears there, a feature that makes the numerical evaluation of the integral easier. The advantage of using identity (5) was pointed out in [33]. Using relation (5) and measuring time and velocity in units of T and U, the dimensionless Maxey-Riley equation valid for any initial condition becomes where n is a unit vector pointing upwards (against gravity). Here three dimensionless parameters appear, the density parameter where ρ p denotes the particle density, the size parameter S = 2 9 a ratio of the particle's viscous relaxation time and the characteristic time T of the flow, and the dimensionless settling velocity This is negative (positive) for heavy (light) particles, called aerosols (bubbles). We shall see that the size parameter S is a more appropriate number to characterize memory effects than the traditional Stokes number, which is St = S/R in our notation. An important condition for the validity of equation (6) is that the particle Reynolds number remains small during the entire dynamics [16]. In addition, S must be small (i.e. the particle's characteristic time scale is much smaller than that of the flow) and the particle radius must be much smaller than the characteristic linear scale L of the flow: a ≪ L.
The last condition assures that the so-called Faxén corrections are negligible [16]. In a recent manuscript Farazmand and Haller [38] prove the local existence of the solutions to (6). With initial particle velocity v(t 0 ) matching the fluid velocity, the solution is twice continuously differentiable in time, but otherwise it might be only differentiable for t > t 0 .
Several attempts have been made to extend the Maxey-Riley equation to finite particle Reynolds numbers (up to a few hundreds) by modifying the particular form of the forces [39,40,41]. Part of all of these approaches is a different form of the history force and in some cases also an empirical non-linear expression for the drag (see [42] for a review). Concerning the range of applicability of the Maxey-Riley equation, it is interesting to note that Maxey and Wang [43] found that up to Re p ≈ 17 equation (1) "may be quite adequate in practice even though not justified by theory". In an experimental study [23] considering particle Reynolds numbers up to 0.5, the standard form of the history force has been found to remain valid. In our simulations the average particle Reynolds number is on the order of unity (see sec. 10). Therefore, we expect the standard Maxey-Riley equation to be adequate for the advection of inertial particles in our case and decide not to consider any finite Reynolds number corrections to the Maxey-Riley equation. In this study we focus on the changes induced by the history force in the motion of inertial particles. Therefore we do not consider effects like e.g., the lift-force.
Equation (6) is a second-order integro-differential equation for the trajectory r(t) of an inertial particle. To any initial condition (r 0 , v 0 ) at time t 0 there is a unique trajectory, but the transition between infinitesimally close neighboring time instants t and t + dt does not only depend on the state q(t) ≡ (r(t), v(t)), but also on all the previous instants, since the kernel of the integral decays rather slowly. Considering a stroboscopic map taken at integer multiples of some time unit, we obtain a map whose value at time n + 1 depends on its value at all previous integers back to the initial time instant: From a dynamical systems point of view, the dynamics is thus that of a non-autonomous map with memory whose range is increasing as time goes on. This is a non-standard problem in which novel features can show up. It is of particular interest which kind of chaos and attractors can be present in such systems.

Estimating the relevance of the history force
To estimate when the history force is important, we compare its magnitude to the Stokes drag (both are viscous effects). From the dimensionless Maxey-Riley equation (6), the ratio of these two forces is seen to scale with √ S. Denoting the typical magnitude of the history force and the drag by F h and F d respectively, we can write where α represents the ratio The value of α might depend on the process in question and is expected to be about unity (see also the discussion in Section 15). Eq. (12) provides an estimate for the magnitude of the history force relative to the Stokes drag. For practical purposes we might consider the history force to be negligible if its magnitude is smaller then 1% of the Stokes drag, i.e. if S < 10 −4 (2/9)/α 2 . It might appear that 1% is already quite small. However, note that we are making order of magnitude estimates here and could be easily off by a factor of 10. In terms of the particle radius a in (8) this 1%-condition means In order to relate this condition to the Reynolds number Re = UL/ν of the flow, we define the Strouhal number Sl as the ratio of L/U and the characteristic time T of the flow Sl = L UT (15) and find a L < 1 as a condition under which memory effects can be neglected. Note that the above argument assumes a smooth flow. For multi-scale flows (turbulence) the length scale L and the parameters Re and Sl in (16) should be characteristic of the small scales (see also the discussion in Section 15). For example, in a smooth flow with Re = 100 and Sl ≈ 1 the particle size has to be about 1000α times smaller then the characteristic length scale of the flow to be able to neglect the history force with some confidence. For 100 µm particles in a case with α ≈ 1 this is not fulfilled in small-scale flows with L < 10 cm.
Another approach is to consider the weakly inertial limit, where S is assumed to be a small expansion parameter. Following Druzhinin and Ostrovsky [21] we obtain the following expression for the velocity of a nearly ideal particle: This is a first-order equation for the trajectory and we see that there are two correction terms to the passive tracer dynamics: the fluid acceleration and the history of the fluid acceleration. Both corrections are proportional to S/R. This ratio determines thus the importance of inertial corrections. Provided that Du/Dt and the derivative of the integral within the parenthesis are of the same order, the importance of memory is determined by √ S, as also found in the previous consideration (when weak inertia was not assumed). Therefore, estimates (14), (16) also hold in the weakly inertial limit (when, of course, the definition of α should be adjusted according to (17)).
Note that the above considerations do not depend on the density of the particle. This is because the density parameter R can be brought to the left-hand side of the Maxey-Riley equation, as has been done in (6). It, therefore, does not influence (at a fixed value of W ) the relative importance of the forces acting on a particle of fixed size.
Thus, the weight of the history force is set by the size parameter S, independent of the particle's density. The density parameter R, however, does influence the importance of inertial effects, as S/R is the prefactor of the correction in (17). The natural parameter to estimate the importance of inertial effects is thus the Stokes number S/R ≡ St, and the size parameter S is the natural quantity to estimate the importance of the history force.
We would like to say a few words about the case of very heavy particles, e.g. water droplets in air, R ≈ 10 −3 . All our previous considerations apply to this case as well.
A new feature of heavy aerosols is, however, that inertial effects might be important in spite of the smallness of S since R is also small. In such cases the pressure term 3R/2 Du/Dt in (6) can always be neglected, but not the others. The relative weight of the history force remains given by α √ S 9/2 and since S is small, memory effects are not expected to be very important. For S ∼ R ≈ 10 −3 , α ≈ 1, for example, the effect is estimated to be 7% of the Stokes drag. When S < 10 −4 , is an (at least) 2%-accurate advection equation for heavy aerosols, often used in the literature, see e.g. [44,14,45]. For even smaller heavy particles: S < < R, the equation simplifies to the passive advection equation with settling: v = u + W n.
Finally, we mention that the requirement of small particle Reynolds number (10) provides another inequality for the particle size. For memory effects to be relevant, inequality (16) should hold with a reversed sign and this sets a lower bound to a/L. The conditions that Re p should not become too large provides an upper bound on the dimensionless particle size.

Numerical scheme
The history force poses the main difficulty for a numerical integration of (6). The first major problem is the singularity of the kernel 1/ √ t − τ at the upper end of the integration, at τ = t. Although the singularity is integrable and the history force is well defined, the numerical evaluation of the kernel near it leads to large errors when standard quadrature schemes are used. Recently, a custom quadrature scheme has been developed [33] where the kernel is treated analytically. It is expressed as a weighted sumˆt where coefficients µ i contain the effect of the kernel. The derivation and specification of the coefficients for higher orders is quite involved [33]. This quadrature scheme is then incorporated in a linear multistep method of the Adams-Bashforth type, where relation (5) is used to evaluate an integral of the history force. The quadrature scheme for the history integral, as well as the whole integration scheme for the Maxey-Riley equation, has been tested in cases where analytical solutions are available, see [33], and the order of convergence has been measured and verified in these cases. The numerical results shown here are obtained with a third-order scheme, i.e. the one-step error is proportional to (∆t) 4 where ∆t is the time step.
The second major problem is the high computational costs for a numerical integration, stemming from the necessity to recompute the history force -an integral over all previous times steps -for every new time step. Therefore the computational costs grow with the square of the number of time steps and can become quite substantial for long integration periods. For the integration of particle ensembles we have addressed this difficulty by parallel computation (on ca. 100 CPUs).

Model flow
An accurate evaluation of the history term is computationally rather intensive, and we would like to use particle ensembles of about 10 6 particles in order to have reliable statistics. Therefore, we choose an analytically given model flow to limit the numerical efforts. As a paradigmatic example, we take the kinematic model of a rather typical fluid phenomenon, the time-periodic shedding of von Kármán vortices in the wake of a cylinder as described in [46] (sometimes called the JTZ flow). This analytical twodimensional flow was shown to faithfully represent the Navier-Stokes dynamics at a (radius-based) Reynolds number Re ≈ 125 and has since then been used in several studies of passive [47,48,49,50,51,52] and inertial [12,53,54,55] chaotic advection in environmental flows. In the JTZ model two vortices are present in the wake at any time (a generalization to more than two vortices has been given recently in [56]). The flow is from left to right.
The period T of vortex shedding is traditionally chosen as the time unit, and the radius of the cylinder as the length unit L. The velocity of the vortex cores is on the order of L/T , whereas the inflow velocity is a factor 14 times larger and is taken as the velocity unit U. Thus the Strouhal number (15) is Sl = 1/14 . One important further parameter is the dimensionless strength of the vortices w = 192/π (for comparison with The Strouhal number in the von Kármán problem is often used in its diameter-based form which corresponds then to 1/7. [53] we shall also use the value w = 24). The stream-function of the model is taken from the literature [53]. Since the flow is two-dimensional, we do not consider here the effect of gravity and hence W = 0 in (6). For a recent investigation of the effect of the history force on sedimentation of aerosols and on rising of bubbles we refer to [57].
Throughout the simulations presented in this paper the initial velocity difference between particle and fluid is zero, i.e. v(t 0 ) − u(t 0 , r 0 ) = 0. Furthermore, most simulations are performed with the starting time t 0 = 0.2 (time zero corresponds to the instant when a vortex is born along the upper right edge of the cylinder). This value has been chosen to minimize the collisions between particles and the cylinder. As the fraction of colliding particles is reasonably small, we have avoided the inclusion of an extra collision rule on the surface and decided to stop integration in such cases and considered the particle to be escaped.
When studying the dynamics of inertial particles without memory, a stroboscopic map can be defined, taken at time instants which are integer multiples of the period T of the flow. This is an autonomous map, in which periodic orbits and invariant sets are defined as usual. In the presence of the memory term, the stroboscopic map becomes of the form of (11), in which periodicity after a finite amount of time is hardly possible, since the state of any time instant is influenced by the entire history. It is therefore a kind of surprise that periodic attractors might exist, as will be seen later.

Choice of particle parameters
The density parameter (7) ranges from 0 to 2. Throughout the paper we shall use two illustrative values R = 0.4 and R = 2 typical for aerosols and bubbles, respectively. The former corresponds to a density ratio ρ p /ρ f = 2 and can be considered to represent sand grains in water. Since air is about thousand times lighter than water, the value of R = 2 can be seen to correspond to the case of spherical air bubbles in water.
For completeness, we mention that water droplets in air are described by a density parameter R = 1.2 × 10 −3 . This is an example of the case of heavy aerosols discussed briefly in Sect. 3. For oil droplets (bubbles) in water we obtain with the density ρ p = 0.9 g/cm 3 of oil R = 0.7. This is close to R = 2/3 of neutrally buoyant particles.
Next we point out that the instantaneous particle Reynolds number (10) can be determined from the data of our model flow. Replacing a and ν in (10) by using (8) and ν = LU/Re we obtain where |v − u| is the dimensionless slip velocity (measured in units of U). This will be used to check the validity of the basic equation (6). The value of the size parameter S will be changed in the range [0.005, 0.09], whereas Re and Sl are set by the flow (to 125 and 1/14, respectively).  Figure 1. Trajectory of a single heavy particle (R = 0.4 -sand in water, S = 0.01) with respect to different dynamics: with memory (red), without memory (blue), and as an ideal tracer (green); initiated at t 0 = 0.2, r 0 = (−1, −1.2) with vanishing slip velocity. The forces acting on the particle in the presence of memory are shown as arrows of different colors (red being the history force). The pressure term is presented with only a third of its magnitude.

Sample trajectories
First, we consider trajectories starting with the same initial condition, and compare different dynamics. Figure 1 shows the inertial trajectory of an aerosol with and without memory along with the trajectory of an ideal tracer (which follows the fluid velocity exactly). One immediately observes that the trajectories differ and that (at the beginning) the trajectory with memory is in-between the other two. The acceleration of the particle can be decomposed into the three contributions on the right-hand side of (6). These terms, called the pressure term, the Stokes drag and the history force are plotted at subsequent time instants in Fig. 1. As the caption indicates, one third of the pressure term is given, for better visibility. We thus clearly see that the pressure term always dominates the other two, but the history force is of the same order as the Stokes drag. A comparison of the full trajectory with that obtained without the memory term (red and blue lines in Fig. 1, respectively) shows that the direction of separation (or approach) of theses trajectories coincides with the direction of the history force. For bubbles, similar findings have been reported in [34].
We also carried out simulations with equation of motion (17) valid in the weakly inertial limit. The results (not shown in Fig. 1) indicate a strong deviation from the true trajectory (red line), as strong as the deviation between this curve and any other curves in Fig. 1. This shows that a value as small as S = 0.01 of the size parameter is not yet small enough to ensure the weakly inertial limit. At the same time, this is a sign indicating that dynamics (17) does not define a slow manifold [54,55] for this S (although it is expected to do so for sufficiently small S values).

Ensemble and escape dynamics
Next, we turn to the dynamics of particle ensembles. A large number N 0 ≈ 1.8 · 10 6 of particles is distributed uniformly in the domain [−1, 4] × [−2, 2] around the cylinder at t 0 = 0.2. All particles are followed up to a certain time when we plot their position. The results for aerosols are shown in Fig. 2. The two distributions differ both in largeand small-scale structures. A characteristic feature of the case without the history force are caustics [4], i.e., the intersection of different branches in the configuration space. This is due to the fact that what we see is a projection of the full pattern in the fourdimensional (x, y, v x, v y ) phase space to a plane. Figure 2 shows that the frequency of occurrence of caustics seems to be lower with memory. We have found similar results for other values of R: 0.3, 0.5, 1.0, 1.2, 1.5, 1.7, 2.0. This effect of the history force is even stronger for bubbles [34]. A quantity closely related to caustics is the collision rate. In accordance to the suppression of caustics, the history force has been found in [34] to lead to a reduction of the collision rate.
In Fig. 2 a filamentary pattern can be seen both with and without the history force. This is in itself an indication for the chaoticity of the advection dynamics in both cases. Since the problem is dissipative, chaotic attractors might be present. Irrespective of their existence, chaos is always present in the form of transient chaos with an underlying chaotic set, a chaotic saddle [58,59]. In the transient case all particles eventually escape downstream and the decay dynamics is best followed by monitoring the total number N(t) of particles not yet escaped a given region up to time t. We distribute N 0 ≈ 1.5 · 10 6 particles uniformly in the domain [0.6, 4] × [−2, 2] outside the cylinder at initial time t 0 = 0.2. A particle is considered to have escaped if it crosses the line x = 5 or it enters a circle of radius r = 1.014 around the cylinder. The using of a circle larger in radius than the cylinder is motivated by excluding very slow (non-hyperbolic) decay characterizing the boundary layer, as done in [12,53]. This definition of "escape" will be used throughout the paper, unless stated otherwise (as, e.g., in Sec. 9). The number N(t) of survivals is then determined as a function of time t. We find monotonous exponential decays with different characteristics for different cases.
When the particles are heavier than the fluid (aerosols) typically all of them eventually escape; an example is shown in Fig. 3a. We see that the emptying process is slower with memory for aerosols. For bubbles attractors can appear for certain parameter values. An attractor presents itself as a plateau in N(t), i.e. a fraction of all the particles never leaves the wake of the cylinder. See for example the case S = 0.04 without memory in Fig. 3b; with memory the plateau, and thus the attractor, disappears. Here we see a profound change in the dynamics of inertial particles induced by the history force: the disappearance of attractors. When there are no attractors in the bubble case, the history force leads to a faster emptying of the wake. In the case of S = 0.02 of Fig. 3b the time for a complete emptying roughly halves when memory is included. Note that the influence of the history force on the speed of the emptying process is opposite for aerosols and bubbles: memory slows down the emptying for aerosols and speeds it up for bubbles.
A quantitative measure characterizing the escape dynamics from the wake is the escape rate [58,59]. It is the rate κ of the exponential decay which sets in after an initial transient time, see Fig. 3. In the case of a plateau in N(t) (due to the presence of an attractor) we set the escape rate formally to zero. Fig. 4 shows the escape rate with and without memory as a function of the parameters R and S. The red line marks the case of a neutrally buoyant particles (R = 2/3), where the escape rate does not depend on S; indeed it is equal to the one of ideal tracers (κ ideal = 0.364). We did not find any influence of memory on the behavior of neutrally buoyant particles. They behave as ideal tracers in this flow when the initial slip velocity is zero. Attractors are quite ubiquitous without memory, see Fig. 4b. One even finds attractors in the aerosol range R < 2/3. Since heavy particles are pushed outside of closed orbits due to the centrifugal force, and intend therefore to escape, attractors are expected to be atypical for aerosols. In special cases they were pointed out to exist [60] and our results provide further examples of such cases. A striking effect of the history force is that it kills attractors at the parameters where they exist without memory. Only very few attractors (none for aerosols) are found with memory.
The escape rates are typically closer to that of the passive (neurally buoyant) case with memory than without. In other words, the speed of the emptying process becomes less R-dependent. Figure 5 illustrates this tendency along two cuts (at R = 0.4 and R = 2). It is clear that the escape rates of aerosols (bubbles) become smaller (larger) due to memory, and that they are typically above (below) that of the passive case. There might be, however, exceptions as the bubble case with S ≤ 0.01 illustrates where the escape rates slightly exceeds κ ideal . Note that the S-dependence is also weaker with memory. The escape process depends, thus, less on both parameters R and S when memory is taken into account.
For completeness, we also present the escape rates on the plane of the Stokes number St = S/R and the density ratio R, Fig. 6. Since in this representation the density plays a role along both axes, the graph of the escape rate is different but the general content remains the same. An important effect of the history force can be observed when studying particles moving along general curved trajectories: the history force is found to counteract the centrifugal force. Heavy particles are pushed outwards of vortex centers by the centrifugal force. Memory is seen to generate a dissipative effect which provides a force pointing towards the center of a vortex, as illustrated by 7a. For bubbles the effect is opposite, i.e. the history force points outwards the center of a vortex; it counteracts the (anti)centrifugal force which pushes bubbles inside of vortices. This behavior has been also found in [61] for the analytically treatable case of a single vortex with the flow field u(r) = |r| e ϕ . The effect visualized in Fig. 7 provides a qualitative explanation for the decrease (increase) of the escape rate for heavy (light) inertial particles in the presence for the bubble. The scale of the arrows for the bubble case is 10 time larger, i.e. in this case the history force is much stronger (this is also true for the the other forces). Note that the history force counteracts the centrifugal force in both cases. of memory. We can summarize that, in general, the history force tends to weaken the difference between particle and fluid and "pushes" the dynamics towards that of the ideal case.

Residence time and nonhyperbolic behavior
A further characterization of the escape dynamics is via the residence time distribution. To any initial position in the flow, we determine the time the particle takes to escape and indicate this time via color-coding. The trajectories are integrated up to a maximum of 100 time units, thus the largest residence time shown is 100. In this section we do not cut out the circle of radius 1.014 along the cylinder surface (unlike in the previous sections), and the escape condition is simply x > 5, in order to also see the effect of the stickiness of the wall. Fig. 8 shows the result with and without memory for bubbles. Very long lifetimes, longer than 100 time units, would indicate basins of attraction of an attractor. In this case no attractors exist and lifetimes on the order of 100 mark slowly moving trajectories which will spend a long time in the boundary layer close to the cylinder surface. We see clearly that the area of initial conditions with long lifetimes (redish colors) is much lower with memory than without. Thus the surface of the cylinder seems to be less sticky when memory is included. Furthermore the structure of the residence time distribution shows a clear dependence on the presence of memory also away from the cylinder surface, e.g. in the region [1.5, 3.5] × [−1.5, 0] where a vortical pattern is present with memory.
In dynamical system terms, the sticky surface of the cylinder leads to a nonhyperbolic behavior marked by a long-term power law decay of the number of survivors. For trajectory ensembles containing initial points close to the surface one finds that the number N(t) of non-escaped particles up to time t behaves for t ≫ 1 as where σ is the exponent of the non-hyperbolic decay [62]. The number of survivors in the range of 100 and 1000 time units exhibits clear straight lines on a log-log plot in both cases, as can be seen in Fig. 9. The slopes are, however, different. We have found a change from σ = 0.55 to σ = 1.21 when memory is included. (For ideal tracers σ = 2 [46].) For bubbles, memory effects lead thus to a faster decay (to an increased decay exponent), weaken non-hyperbolicity, but do not destroy the power-law behavior. For aerosols, the effect of the cylinder wall appears to be much weaker both with and without memory, and we don't find any change of the non-hyperbolic decay exponent.

Statistics of forces, slip velocity and acceleration
In section 7 we briefly discussed the forces acting along a sample trajectory. To be able to make more general statements we now turn to the statistics of these forces. we sampled over all time instants at which a particle has not yet escaped the wake of the cylinder. Fig. 10a shows the probability density function (PDF) of the magnitude of the three forces appearing on the right-hand side of (6). Let us first compare the averages of the different forces, shown as vertical lines in Fig. 10a. We see that our observations from Sec. 7 based on a single example trajectory also hold true for the averages: the pressure term is the dominant one and the history force (red curve in Fig. 10a) is of similar magnitude as the Stokes drag. Looking at the tails of the PDFs we observe that the pressure term has a particularly large amount of extreme events and the history force has a slightly longer tail (more extreme events) then the Stokes drag.
The ratio of the averages of the history force and of the Stokes drag in Fig. 10a is F h /F d = 1.21. In view of this finding, (12), and the fact that √ S = 0.1, we conclude that parameter α is 2/9 × 12.1 = 5.7 in this setting. The observation that the α value is not close to 1 shows that the importance of the history force cannot be estimated solely by √ S; parameter α is also relevant. If we neglected α, we would underestimate the history force by a considerable amount.
In order to be a useful parameter for estimating the importance of the history force, α should be independent of S. Indeed, for R = 0.4 we find α ∈ [4.7, 6.0] and for R = 2.0 α ∈ [7.2, 7.7] when S is varied between 0.05 and 0.005. For the purpose of rough estimates we can thus consider α constant, with a value around 5 (in this flow). Fig. 10b shows the PDF of the magnitude of the slip velocity |v − u|. It is important to see that the average is much below 1 in both cases, indicating that the deviation between the particle and the flow velocity is typically much smaller than the flow velocity itself. Furthermore, the average dimensionless slip velocity 0.153 of the memoryless case is reduced to 0.085 (to ca. 55%) when the history force is included. Thus, particles with memory follow the fluid more closely. In harmony with this, the tail of the slip velocity PDF becomes shorter with memory, i.e. there are less strong "detachments" from the fluid flow. We find the same tendency for bubbles albeit somewhat weaker, e.g. for R = 2, S = 0.01 the slip velocity is reduced to 85% by memory. All this illustrates that the history force provides an important viscous contribution so that the resultant of it and the drag ensures a faster relaxation towards the flow velocity than the drag alone.
Equation (20) shows that the slip-velocity is proportional to the instantaneous Re p . Accordingly, the top axis of Fig. 10b displays Re p . Its average for particles with (without) memory is found to be 0.76 (1.36). The assumption that the particle Reynolds number remains of order unity or smaller is thus found to be fulfilled on average, with a somewhat sharper bound with memory.
The statistics of the acceleration of particles (not shown) is also consistent with the picture above. Ideal tracers have the highest average acceleration and the most extreme events (the longest tail in the PDF). The case without memory exhibits the lowest accelerations, and the inclusion of memory increases the average acceleration, as well as, the probability of extreme events. Thus we see again that the statistics come closer to the ideal case when memory is included.

Periodic attractors
Let us now turn to a parameter set where attractors are present both with and without memory: w = 24 (instead of 192/π), where w describes the strength of the vortices behind the cylinder, and R = 1.5, S = 0.03. The residence times are depicted in Fig. 11. Their structure changes significantly when memory is taken into account, as we have already seen in section 9. Here we included again the circle of radius 1.014 as an escape condition to exclude the non-hyperbolic decay. Thus residence times of 100 (the maximum integration time) show initial condition leading to an attractor. In the presence of memory, the basin of attraction (marked by red) becomes significantly smaller; its area is reduced by a factor of about 2.5. Thus, we can say that memory makes the attractor less attractive in this case. We note, that for this choice of w there are also trajectories of ideal tracers which never leave the wake. Tiny integrable regions exist which appear as red dots in Fig. 11c. Comparing the residence time plots we see again that the case with memory comes closer to that of ideal tracers.
It is worth comparing the attractors with and without memory. Both attracting objects appear to be periodic orbits of period-1 (i.e., with the same period as the vortex shedding period), and are of similar form, as Fig. 12a shows. Note that the history force does not vanish on the attractor with memory. A basic difference between the two cases is that the attractor of the memoryless case is stationary by the time the plot is made, but the cycle-looking object seen in the presence of memory keeps changing very slowly in time. It is thus a 'quasi-attractor' only. A real attractor is expected to be reached (even in a numerical sense) after a very long time only.
We examine the convergence of the trajectory to the attractor through the quantity |x(n) − x * |, where x(n) and x * is the x-coordinate of an approaching trajectory after n time units and of the asymptotic attractor (x * = lim n→∞ x(n)), respectively. This corresponds to taking snapshots at integer multiples n of the period, i.e. a stroboscopic map. The time dependence of this quantity is shown in Fig. 12b for the case with memory. The convergence follows a t −1/2 law, where the exponent is obviously related to the power law behavior of the kernel of the history force (5). The algebraic convergence is in strong contrast to normal dynamical systems where the convergence is exponential. This means that with memory the convergence is very slow and that there is no characteristic time scale characterizing this convergence.
The fact that a periodic attractor is possible with memory can be explained by the observation that after a very long time the memory of the approach to the attractor decays away, and the dynamics remembers only what has been in the close vicinity of the attractor, on a loop like the continuous line in Fig. 12a, and a convergence becomes thus possible.

Interpretation in terms of snapshot attractors
Using this classical single trajectory picture, one concludes that the attractor can be reached after a very long period of time only. There is, however, an alternative view, that of particle ensembles, also available. In autonomous systems the two views are equivalent, which is not so obvious in non-autonomous problems, like ours (see (11)). In this class, one can then define a snapshot attractor as an object which attracts all the trajectories initialized in the infinitely remote past within a basin of attraction [63,64]. It can be obtained by monitoring an ensemble of particle trajectories all subject to the same non-autonomous equation of motion. After a characteristic dissipative time, the ensemble traces out a snapshot attractor. This attractor might, however, move continuously in time.
In the dynamical systems community, the concept has been known for many years [63]. The idea proved to be particularly well suited for understanding the advection of passive particles in random flows [65,66,67], and explains experimental findings [68]. More recently, snapshot attractors have turned out to be promising tools to describe the variability of climate dynamics in a novel way [64,69,70,71,72]. In these settings the non-autonomous driving is typically a random noise or a sustained chaotic process. The snapshot attractor is then ever changing, and appears to be a fractal object whose shape is evolving in time. In cases with an eventually vanishing driving, the snapshot attractor might have a time-dependence that ceases asymptotically.
For our problem of chaotic advection with memory, where a slow convergence towards a traditional attractor takes place, one can assume that close to this attractor the effect of the history force for neighboring members of the particle ensemble can locally be considered as a non-autonomous perturbation superimposed on the usual inertial particle dynamics. An ensemble of particles might then converge to a snapshot attractor, a fixed point on the stroboscopic map, after some time. This snapshot fixed point attractor is then slowly shifted towards the asymptotic traditional attractor. Here a separation of time scales is expected to occur since the convergence to the snapshot attractor is a usual dissipative effect, and hence this decay should be exponential, whereas the slow shift is due mainly to the diffusive kernel, and follows thus a power law.
To verify this we choose an ensemble of 100 particles approaching the periodic attractor, and compute the center of mass r(n) (i.e., the average of the spatial location vector over the ensemble) and the radius of the ensemble |r(n) − r(n) | at integer times n. The center of mass represents the time-dependent position of the snapshot attractor whereas the radius is a measure of convergence to the snapshot attractor. The center of mass is found to exhibit the same power-law behavior, of power −1/2, as x(n). However, we observe that the radius of the ensemble has an initial roughly exponential decay, see Fig. 13a, which crosses over to an algebraic decay of power −3/2. Thus we conclude that the point-like object formed after about 70 time units is a snapshot fixed point attractor emerging as an effect of memory.
To deepen this picture, we plot in Fig. 14 the center of mass and the radius of the ensemble in the plane of the fluid at integer times. Fig. 14a shows time instants between n = 200 and 300, whereas Fig. 14b covers the range (200, 1000). The scale is strongly magnified but the plot clearly indicates that the fixpoint snapshot attractor slowly moves. Even after 1000 time units it is still away from the limiting location. The latter is determined from a self-consistent fit to the algebraic decay process of the center of mass.
The evolution can thus be split into a short-time and a long-time regime. In the first one, up to about 70 time units, a convergence to a snapshot attractor takes place. In the second one (for n > 70) a snapshot attractor is reached but it is shifted slowly towards an asymptotic fixed point, the traditional attractor.

Chaotic attractors
Chaotic attractors with memory are very rare in the parameter range investigated and we had to perform an extended parameter-scan to find one. In Fig. 15a we see an ensemble of 100 trajectories (in continuous time) tracing out this chaotic attractor with memory. A magnification, Fig. 15b, indicates clear fractal properties. The Lyapunov exponent of this attractor is λ = 0.18 ± 0.01. It is positive, which shows that the attractor is chaotic indeed. Details of how the largest Lyapunov exponent can be determined in such a system with long memory is given in Appendix A. We note that at this parameter values (given in the caption to Fig. 15) the attractor without memory is not chaotic. Fig. 15c exhibits the stroboscopic picture of the chaotic attractor with memory. It shows a magnification in two different time intervals: t ∈ [600, 800] and t ∈ [800, 1000]. A shift on the order of 10 −4 units can be observed, a clear manifestation of the fact that this chaotic attractor is slowly drifting, similar to the fixed point snapshot attractor discussed in the previous section.

Chaotic saddles
In contrast to permanent chaos, transient chaos is ubiquitous in the problem as the typical form of chaos in this open flow. Chaotic transients are found at any parameter, irrespective of the existence of attractors. The degree of chaos can thus best be compared by comparing the underlying chaotic saddles which exist both for ideal tracers and inertial particles with or without memory.
To construct the saddle, we consider particles with long residence times, longer than 15 time units. The idea is that these particles come close to the saddle and then leave the wake. The longer the particle stays in the wake, the closer it comes to the saddle. One might suppose that at half of their residence time (t = 7) the particles are closest to the saddle and thus their positions at this time trace out the saddle. This is supported by the observation that this set of points at t = 7 is only marginally different from those at t = 6 and t = 8, i.e. the set is (approximately) invariant. Fig. 16a shows the saddle constructed by this method for three different advective dynamics. The initial positions of particles which stay long in the flow and thus come close to the saddle trace out the stable manifold of the chaotic saddle, and are plotted in Fig. 16b.
We see that the saddle and its stable manifold change considerably when the history force is taken into account. It is interesting to note that the saddle and its stable manifold differ much more from those of ideal tracers than from the memoryless case. This illustrates that the general rule according to which the presence of memory makes the dynamics to be similar to that of ideal tracers is not true in all possible aspects.
To further characterize the saddle and the corresponding dynamics we calculated the Lyapunov exponent, as described in Appendix A. The Lyapunov exponent with the history force is found to be λ = 0.91 ± 0.02. This is larger then without memory, when λ = 0.79 ± 0.02, and is close to that of ideal tracers for which λ = 0.92 ± 0.02. Thus particles with memory are dynamically more unstable, and their measure of instability is closer to that of ideal tracers.

Summary and discussion
In summary, memory effects can have a very significant influence on inertial particles. In the presence of the history force, compared to the memoryless inertial case, we find: • a rare appearance of attractors, especially of chaotic ones, • a less frequent appearance of caustics, • a decrease of the centrifugal effect, • a slower (faster) escape for aerosols (bubbles), • a weaker effect of the cylinder wall, • a stronger dynamical instability on chaotic sets, • a very slow, algebraic convergence to attractors (when defined in the usual single particle picture), • a tendency of the particles to follow the flow more closely.
As a trend, the dynamics become more similar to that of ideal tracers when the history force is included (although exceptions can also be found, see, e.g., Fig. 16).
In one principal aspect the dynamics basically differs, however, from both the usual inertial and the passive tracer dynamics. Due to memory, the problem is strictly speaking high-dimensional, and methods known from low-dimensional dynamical systems are not necessarily applicable. For understanding the slow convergence towards periodic attractors, the concept of snapshot attractors appears to be useful. More generally, we assume that the snapshot interpretation, i.e., considering the instantaneous position of an ensemble of trajectories started in the remote past, is applicable for any kind of asymptotic sets (chaotic attractors, chaotic saddles) in advection with memory effects. In this new picture, an exponential convergence to a snapshot-like attractor is expected to show up, which is then slowly changing afterwards.
Let us come back to estimation (12) of the importance of the history force. A question of relevance is whether one can determine a typical value for α in general. For the following argument we assume the time to be dimensional again. Imagine a signal f (t) which varies on the time scale τ . Then we can estimate the magnitude of the derivative df dt by |f |/τ . Extrapolating this reasoning to fractional derivatives, we expect the magnitude of d dt 1/2 f to be roughly 1/τ |f |. Thus, setting f = v − u in (13) we In our flow the time unit T is chosen, traditionally, as the shedding period of the vortices.
As particles "live" on the small scales, τ should be the smallest time scale of the flow + .
In the von Kármán problem there are two characteristic times: the shedding period T and the convective time scale L/U. The later is the smaller one and can be considered to be the time a particle needs for going around a vortex. We thus take τ = L/U in (23) and find, in view of (15), that Using the value Sl = 1/14 we obtain α ≈ 3.7, which is quite close to the measured range 4.7 ≤ α ≤ 7.7.
The reason α turns out to depend on Sl is that we have used the shedding period T to nondimensionalize the equation of motion, rather than the smaller time scale τ = L/U. Thus we conclude, a posteriori, that a more natural choice is to use the convective time scale from the very beginning. Denoting the corresponding quantities in this time unit with a prime, we obtain α ′ ≈ 1 and with a modified size parameter S ′ = 2 9 a 2 /ν τ .
(26) ¶ The factor √ T appears because time is dimensional here in contrast to (13). + In principle the typical time scale of v − u depends on the particle parameters. Nevertheless, it should be close to the small time scale of the flow. We choose this time scale for our order of magnitude estimates. This choice is further supported by the following finding: using (17) one can obtain a Taylor-expansion of α (13) in S, where the zeroth order approximation α ≈ D 1/2 Dt is solely a property of the flow (here (D/Dt) 1/2 is defined in the same way as (d/dt) 1/2 but with the derivative taken along a fluid element trajectory).
We thus see that for the size parameter (and the Stokes number) -a Lagrangian quantity -the use of the smallest time scale, the convective time scale L/U, is natural. It leads to simpler expressions, in contrast to using the shedding period T which is the usual choice in the Eulerian characterization of this flow.
Let us come to a more general statement about other flows. From (25) and (26) we obtain the rather simple formula for the estimation of the history force relative to the Stokes drag. We expect this to be a good estimation also in other flows as long as the time scale τ is chosen to be the smallest time scale of the flow. This choice becomes particularly important in multiscale flows. Note also that the convective time scale need not necessarily be the smallest time scale in other flows. However, when this is fulfilled and when the flow is smooth with only one characteristic length scale L (like the flow examined in this paper) we get Here it is interesting to note that the Faxén corrections to the Stokes drag go as (a/L) 2 [16]. The corrections due to memory are thus more important than the Faxén corrections when a < L. Let us finally remark on the possible role of memory effects in turbulent flows. In turbulence, the time scale τ should be the Kolmogorov time scale τ K [73]. Using ντ K = η 2 in (27), where η is the Kolmogorov length, we obtain as an estimate for the importance of the history force. We thus expect memory effects to be negligible for a < 10 −2 η. For a around 0.1η the history force should be around 10% of the Stokes drag and is thus expected to be clearly observable. However, when a is close to η, the history force can be as important as the Stokes drag indicating strong memory effects even in turbulent flows * . (Note that in such cases Faxén corrections might also be important). (iii) repeat step 1 and 2 until the desired integration time t end is reached This rescaling makes sense because q ′ and q are rather close and thus ∆q evolves according to the linear equation (A.1). Therefore, at any iteration of the above procedure, ∆q(t)∆q min /∆q max is an (approximate) solution of (A.1) and q ′ (t) := q(t) + ∆q(t)∆q min /∆q max is a valid perturbed trajectory.

Acknowledgments
The key idea is that when the upper bound ∆q max is reached, we just choose another solution q ′ which starts at a smaller initial perturbation ∆q 0 ; the new perturbation is smaller but proportional to the previous one. As ∆q evolves according to a linear equation, we can efficiently obtain the new solution by simple rescaling. Because we rescale the previous solution and continue with the integration only fort < t, we avoid computations with too small ∆q.
With the history force it is important to adjust q ′ for the whole past t 0 ≤ t ≤t because the evolution depends on the whole past. Without memory one would only need to adjust q(t), which leads to a standard method of computing the Lyapunov exponent [75,58].
When a perturbation trajectory ∆q has been computed up to t = t end , we obtain the Lyapunov exponent from where t end is a large integration time.
In section 14 we are dealing with transient chaos. Therefore the length of the trajectories is limited and it is difficult to accurately compute the Lyapunov exponent from a single trajectory. Therefore we use a preselected ensemble of trajectories which stay close to the saddle for a sufficiently long time, as in other cases of transient chaos [62]. Equation (A.6) defines the Lyapunov exponent as the slope of the function ln |∆q(t)|, computed using its values at t = t 0 and t = t end . We extend this idea to an ensemble by computing the average and making a linear fit to d(t) to obtain the slope which is the Lyapunov exponent λ.
In section 14 we choose trajectories which stay in the wake for at least 15 time units and make a linear fit in the range t ∈ [5,15].