Fabiola A Gerosa
email: fabiola.gerosa@oca.eu

Héloïse Méheut

Jérémie Bec

Clusters of heavy particles in two-dimensional Keplerian turbulence

Protoplanetary disks are gaseous systems in Keplerian rotation around young stars, known to be turbulent. They include a small fraction of dust from which planets form. In the incremental scenario for planet growth, the formation of kilometer-size objects (planetesimals) from pebbles is a major open question. Clustering of particles is necessary for solids to undergo a local gravitational collapse.

To address this question, the dynamic of inertial particles in turbulent flows with Keplerian rotation and shear is studied. Two-dimensional direct numerical simulations are performed to explore systematically two physical parameters: the rotation rate, which depends on the distance to the star, and the particle response time, which relates to their size. Shear is found to drastically affect the characteristics of the turbulent flow destroying cyclones and favoring the survival of anticyclones. Faster Keplerian rotation enhances clustering of particles. For intermediate sizes, particles concentrate in anticyclones. These clusters form in a hierarchical manner and merge together with time. For other parameter values, solids concentrate on fractal sets that get more singular with rotation. The mass distribution of particles is then found to be multifractal with small dimensions at large orders, intriguing for triggering their gravitational collapse. Such results are promising for a precise description and better understanding of planetesimal formation.

I. INTRODUCTION

A stellar system appears after the gravitational collapse of a molecular cloud. While most matter accumulates in the young star, a tiny fraction (less than a few percents) sediments in a protoplanetary disk, which rotates around the star and serves as a stage for the formation of planets, moons, asteroids, and other small bodies. Such disks are mainly composed of gas, but also contain some sub-micrometer dust particles (≈ 1% in mass) that are the building blocks of larger solid bodies. One of the major open questions in the theory of planet formation is the emergence of planetesimals (objects with sizes of the order of kilometers) that are subsequently able to further accrete solids and even capture the surrounding gas. Indeed classical scenarios, in which dust particles grow in a hierarchical, sequential manner [START_REF] Safronov | Evolution of the protoplanetary cloud and formation of the earth and the planets[END_REF], are facing several difficulties.

A first issue relates to the outcomes of collisions between solids [START_REF] Blum | Aspects of laboratory dust aggregation with relevance to the formation of planetesimals[END_REF]. At sizes larger than centimeters, 1 collisions between pebbles are expected to involve too violent impact velocities. They can hardly lead to coagulation, but rather result in bouncing or fragmentation, making hierarchical growth by far too slow and inefficient at such sizes. A second issue comes from the radial drift that is experienced by solid particles [START_REF] Weidenschilling | Aerodynamics of solid bodies in the solar nebula[END_REF]. Protoplanetary disks are indeed expected to have a large-scale radial pressure gradient, so that the gas rotates with a slightly sub-Keplerian speed. Solid particles are insensitive to this pressure gradient and rotate with the exact Keplerian velocity. They are thus dragged by the gas that slows them down, making them drift radially toward the central star. This process occurs on too short time scales to form planetesimals by sequential growth.

Drift and fragmentation make the growth of planetesimals difficult. Nevertheless, it is generally assumed that such effects can be overcome by a direct self-gravitational collapse of dust particles that can bypass hierarchical coalescence. This requires both gravitational interactions to be dominant over the dispersion induced by the Keplerian shear, and velocity differences between dust grains to be sufficiently feeble. This could occur during the settling of solids in the disk's equatorial plane [START_REF] Goldreich | The formation of planetesimals[END_REF]. We focus here on another scenario involving spatial inhomogeneities, whose presence provides favorable conditions. Self-gravity can indeed become dominant in regions where the density of solids is large and their velocity dispersion narrow [START_REF] Gole | Turbulence regulates the rate of planetesimal formation via gravitational collapse[END_REF]. In such a scenario a major role is played by the turbulence of the carrier gas [START_REF] Ishihara | Dust coagulation regulated by turbulent clustering in protoplanetary disks[END_REF][START_REF] Hartlep | Cascade model for planetesimal formation by turbulent clustering[END_REF]. It is even largely admitted that the back reaction of the particles on the gas could itself trigger turbulence. Underlying mechanisms involve the streaming instability [START_REF] Johansen | Protoplanetary disk turbulence driven by the streaming instability: nonlinear saturation and particle concentration[END_REF] or other resonant drag instabilities, in which the decoupling between gas and particles (through settling, drift, buoyancy, or magnetic waves) strongly amplifies tiny fluctuations in the concentration of solids [START_REF] Squire | Resonant drag instabilities in protoplanetary discs: the streaming instability and new, faster growing instabilities[END_REF]. Recent evidences have however been obtained that, on the one hand, an already developed turbulent state impedes the growth of such instabilities [START_REF] Chen | How efficient is the streaming instability in viscous protoplanetary disks?[END_REF][START_REF] Umurhan | Streaming instability in turbulent protoplanetary disks[END_REF] and that, on the other hand, they prohibitively require starting from already large particle sizes [START_REF] Carrera | The streaming instability cannot form planetesimals from millimeter-size grains in pressure bumps[END_REF].

Besides planet formation, the origin and nature of turbulence in protoplanetary disks, and how it couples to particle dynamics, is key to several other astrophysical questions. For instance, disks are known to have a rather short lifetime (of the order of a few millions years). The dissipation of angular momentum by the gas molecular viscosity is orders of magnitude too low to explain this fast accretion onto the central star. Yet turbulence is an excellent candidate to transport matter on shorter timescales, hence providing estimates on the expected turbulent intensity of disks. Transport mechanisms in protoplanetary disks, possible instabilities and the development of a turbulent state are still the topic of intense work (see, e.g., [START_REF] Armitage | Dynamics of protoplanetary disks[END_REF][START_REF] Lesur | Hydro-, magnetohydro-, and dust-gas dynamics of protoplanetary disks[END_REF] for a review). Whether it originates from magneto-rotational effects or purely compressible or incompressible hydrodynamical instabilities, disk turbulence is confined, in rotation, stratified, and thus generally displays fairly strong two-dimensional features, as for instance long-lived vortical structures [START_REF] Bracco | Spotted disks[END_REF].

The presence of such vortices might clearly impact dust accretion and consequently the formation of planetesimals. While heavy inertial particles are usually expelled from rotating regions by Maxey's centrifugal effect [START_REF] Maxey | The gravitational settling of aerosol particles in homogeneous turbulence and random flow fields[END_REF], the presence of shear and rotation in protoplanetary disks actually has an opposite influence. Numerical [START_REF] Tanga | Forming planetesimals in vortices[END_REF][START_REF] Bracco | Particle aggregation in a turbulent keplerian flow[END_REF], analytical [START_REF] Barge | Did planet formation begin inside persistent gaseous vortices?[END_REF][START_REF] Chavanis | Trapping of dust by coherent vortices in the solar nebula[END_REF] and experimental [START_REF] Solomon | Observation of anomalous diffusion and lévy flights in a two-dimensional rotating flow[END_REF] studies indeed give evidence that dust particles cluster in some of the structures of the flow. More specifically, solid particles experience a Coriolis force that can either accelerate their ejection from cyclones (vortices with the same sign as global rotation) or, when strong enough, overcome centrifugal forces in anticyclones. This pushes the particles toward the core of anticyclonic vortices, possibly creating extremely dense point clusters that are excellent candidates for gravitational instabilities. These mechanisms have been known for years, but their study remained essentially qualitative or limited to model flows. Understanding further what is the impact of vortex clustering on planetesimal formation requires more quantitative insights, such as estimates on involved timescales in the presence of turbulent fluctuations, dependence on physical parameters, nature of the associated mass distribution, etc. The present work aims providing insights on such questions.

We present in this paper 2D direct numerical simulations with the idea of shedding some lights on the complex problem of dust dynamics in turbulent protoplanetary disks. We investigate the clustering properties of solid particles in a forced, developed turbulent incompressible flow subject to Keplerian rotation and shear. We follow the "shearing box" approach, which consists in solving locally the Navier-Stokes equation with a mean constant shear and with periodic boundary conditions on a domain that is distorted in the direction of the mean velocity. A number of simulations are performed varying systematically the two main physical parameters of the problem: the particle response time τ p , which measures their inertia and their lag on the gas flow, and the rotation rate Ω, which prescribes the mean shear. We then analyse quantitatively various dynamical and statistical features of solid particles, and in particular their clustering properties when flow structures, typical of a turbulent flow in rotation, are present. Tools borrowed from the study of dynamical systems, such as Lyapunov exponents and fractal dimensions, are used to characterize and quantify particle clusters.

Before focusing on particles, we find that the mean shear has noticeable impacts on the properties of the turbulent flow. Both energy and enstrophy budgets, as well as energy and enstrophy spectra, change substantially with increasing rotation. Strong rotation also leads to pronounced anisotropies in the flow and thus to increased skewness and kurtosis of the vorticity distribution. Shear is found to cause the preferential formation and survival of anticyclonic vortices. Concerning particles, we find that at specific values of the physical parameters, they can form point clusters located inside anticyclones. This strong clustering requires large-enough values of the rotation rate Ω and intermediate values of the particle response time τ p . These clusters are found to form in a hierarchical manner and to merge one with the other when time increases. We propose a simple kinetic model that catches most aspects of their evolution and distribution. Outside this extreme regime, solid particles can nevertheless concentrate on dynamically evolving fractal sets whose dimensions non-trivially depend on the rotation rate. Evidence is obtained that the mass distribution of particles is then multifractal, with a dimension that decreases and saturates at large orders. Such a behavior is key to quantify the probability to get a large local density of solids and to trigger self-gravitating instabilities. In the strong clustering regime, such instabilities would occur during transients.

This paper is organised as follows. We first introduce and characterise in section II the flow used for the numerical simulations and we explain in section III the considered dynamics of particles. In Section IV we show the results of the simulations, highlighting the different behaviors observed for the distribution of inertial particles. In section V and section VI the mass distribution of particles is inspected in the regime of fractal and strong clustering, respectively. Finally, in section VII we conclude the paper summarizing the results of this study and connecting it back to the astrophysical context.

II. THE TWO-DIMENSIONAL SHEARING BOX

Dust and gas dynamics in protoplanetary disks, as in most astrophysical situations, involve gravity and rotation. The balance between gravity and centrifugal forces defines the Keplerian angular velocity, which is given by

Ω(r) = G M /r 3 ( 1
)
where M is the mass of the central star, G is the gravitational constant, and r is the distance to the central star. This velocity profile corresponds to a rotating flow with a strong shear (see Fig. 1).

FIG. 1. Schematic view of a protoplanetary disk and of the rotating shearing box that we here consider.

A. Model and equations of motion

We focus on a small box whose edge is located at a distance r 0 from the star and which is rotating with the reference angular velocity Ω(r 0 ) (see Fig. 1). We then write the dynamics in the co-moving coordinates x = rr 0 and y = r 0 (θ -Ω(r 0 )t), where θ denotes the azimuthal angle. When the reference distance r 0 is much larger that the box size, curvature terms can be neglected. Local variations of the Keplerian angular velocity are shearing our small box. We indeed have Ω(r) -Ω(r 0 ) ≈ (rr 0 ) d dr Ω(r 0 ) = -3 2 Ω(r 0 ) x/r 0 . In the co-rotating frame, gravity and centrifugal terms balance at r = r 0 . Away from r 0 , a tidal force and a Coriolis term must be added. This model, introduced by Hill [START_REF] Hill | Researches in the lunar theory[END_REF] in the late 70's, leads to write the momentum conservation for the gas as:

∂ t v + v • ∇v = - 1 ρ g ∇p + ν∇ 2 v -2Ω × v + 3Ω 2 x e x , (2)
where v is the total gas velocity in the co-moving frame. ρ g and ν are its mass density and kinematic viscosity, respectively. We have here denoted Ω = Ω(r 0 ) = const and Ω = Ω e z . For an incompressible flow, the equation of continuity reads:

∇ • v = 0, (3)
which hence prescribes the pressure p in (2). The use of incompressibility is justified by the fact that, at the small scales that we consider, density waves are less important. Therefore, density can be assumed constant over time. Note that averaging (2) leads to the mean shear

∂ x v y = -3 2 Ω.
We focus here on two-dimensional flows. Several considerations partly legitimize this choice. A first idea comes from the limit of small Rossby numbers Ro = U/(LΩ)

1, where U and L are characteristic velocity and length scales of the flow. The Taylor-Proudman theorem [START_REF] Proudman | On the motion of solids in a liquid possessing vorticity[END_REF][START_REF] Taylor | Motion of solids in fluids when the flow is not irrotational[END_REF] ensures that any motion that occurs on timescales longer than Ω -1 becomes independent of the z coordinate. The flow is then approximately twodimensional and characterized by the well-known Taylor columns [START_REF] Taylor | The motion of a sphere in a rotating liquid[END_REF], that are long columnar eddies oriented parallel to the rotation axis. A second idea comes from the fact that protoplanetary disks are expected to be highly stratified in the vertical direction z. This could be another mechanism by which the gas dynamics becomes almost two-dimensional. In this case, rather than being long and columnar, eddies are now flat with the shape of pancakes. Their formation originates then from internal gravity waves with group velocity perpendicular to the gravitational force [START_REF] Davidson | Turbulence: an introduction for scientists and engineers[END_REF]. Both kinds of structures, columns and pancakes, typically coexist in the presence of both rotation and stratification [START_REF] Cambon | Turbulence and vortex structures in rotating and stratified flows[END_REF]. However, the global dynamics of protoplanetary disks actually involve baroclinic amplifications, thermal transfers, radial stratification, elliptical instabilities, etc., that are all affecting the formation and survival of coherent eddies [START_REF] Barge | Vortices in stratified protoplanetary disks-from baroclinic instability to vortex layers[END_REF]. Studying such instabilities is beyond the scope of our work. We rather focus on two-dimensional solutions of (2) with turbulent fluctuations and vortex dynamics maintained in a developed, statistically steady state thanks to the addition of an external random forcing.

In our settings, the fluid velocity solves the twodimensional incompressible Navier-Stokes equation. The divergence-free turbulent fluctuations u = v + 3 2 Ω x e y are such that the associated vorticity

ω = ∇ T • u ≡ ∂ x u y -∂ y u x solves ∂ t ω + v • ∇ω = ν∇ 2 ω -α ω + f ω . (4)
Note that, because of incompressibility, the Coriolis force does not contribute to the vorticity dynamics. We are here applying a stochastic forcing f ω in order to maintain a developed turbulent state. This forcing is assumed Gaussian, with zero mean, homogeneous, isotropic, white in time, and with spatial correlations concentrated at large scales. As common in two-dimensional turbulence, we suppose that vorticity experiences a linear friction with coefficient α. Such a term prevents kinetic energy from piling up in the flow. It typically originates, either from the friction between different layers in the stratified case, or from the effect of boundaries when the twodimensional effective dynamics corresponds to a vertical average.

Figures 2a, b, andc show snapshots of the gas vorticity field (normalized by τ -1 ω ) for various values of the rotation frequency Ω. One clearly observes, at a qualitative level, that shear sustains anticyclonic vortices (rotating in the opposite direction of Ω, i.e. with a negative sign, in blue). Because the space average of the fluctuating vorticity ω is conserved and equal to 0, positive values (in red) distribute in the background as denser and denser filaments. Another observation is that shear tends to stretch and align anticyclones with the direction of the mean flow. At large values of the rotation, structures typically look like ellipses. We will turn back later to estimate how their aspect ratio depends on the shear rate. The developed turbulent state that we observed in Fig. 2 is characterized by two inviscid quadratic invariants, namely the kinetic energy E = 1 2 u 2 and the enstrophy Z = 1 2 ω 2 , whose budgets satisfy

dE dt = 0 = -2ν Z -2α E -3 2 Ω R + ε I , (5)
dZ dt = 0 = -2ν P -2α Z + η I . (6)
We have here introduced the Reynolds shear stress R = u x u y and the palinstrophy P = (1/2) ∇ω 2 . The quantities ε I and η I are the injection rates of kinetic energy and enstrophy, respectively, and are fully determined by the stochastic forcing. Because of the absence of vortex stretching in 2D, the only source of enstrophy in [START_REF] Ishihara | Dust coagulation regulated by turbulent clustering in protoplanetary disks[END_REF] is forcing. This implies that, contrarily to 3D, energy is not dissipated by viscosity but is transferred to large scales by an inverse cascade (see, e.g., [START_REF] Boffetta | Two-dimensional turbulence[END_REF]). Still, vorticity gradients, and thus palinstrophy, can become very large in 2D, leading to a direct cascade of enstrophy. The two injection rates define a forcing length scale

f = (ε I /η I ) 1/2
. Similarly, by balancing the different dissipative rates in ( 5) and ( 6), one can introduce a friction scale α = (E/Z) 1/2 and a viscous dissipation scale ν = (Z/P) 1/2 . These length scales define two Reynolds numbers: an outer-scale Reynolds number R α = α / f which measures the ratio between inertial and frictional forces, and a viscous Reynolds number R ν = ( f / ν ) 2 , which balances inertial and viscous forces. These two numbers prescribe the extensions of the inverse energy cascade and of the direct enstrophy cascade, respectively. We here focus on the direct cascade of enstrophy and we prescribe α f ν , so that R α is of the order of 1 and R ν 1. As seen from Fig. 2, a mean shear (Ω = 0) develops anisotropies in the flow and the energy budget is affected by the Reynolds stress. The importance of shear has to be measured by non-dimensionalizing it with a characteristic time scale of the flow. Still, the flow timescales are themselves modified by shear, so this choice cannot be made a priori. In our protocol, the only time scale that is prescribed by the simulation setup is the forcing time scale

τ f = ( 2 f /ε I ) 1/3 = η -1/3 I
. The influence of shear on the energy and enstrophy budget is then measured by the non-dimensional shear rate parameter σ = 3 2 Ω τ f . The three illustrating values shown in Fig. 2 are (a) σ = 0, (b) σ = 5.7, and (c) σ = 14.2. In the developed regime attained once σ has been prescribed, the shear needs to be compared to the typical dynamical timescale of the direct cascade, namely

τ ω = (2Z) -1/2 = ω 2 -1/2 .
Direct numerical simulations are performed using a pseudo-spectral solver. To construct periodic solutions that account for the mean flow, we follow [START_REF] Rogallo | Numerical experiments in homogeneous turbulence[END_REF][START_REF] Rogers | The structure of the vorticity field in homogeneous turbulent flows[END_REF][START_REF] Hawley | Local Three-dimensional Magnetohydrodynamic Simulations of Accretion Disks[END_REF][START_REF] Pumir | Turbulence in homogeneous shear flows[END_REF] and integrate (4) on a distorting frame defined by x = x, y = y + 3 2 Ω t x. The integration domain is the periodic rectangle (x , y ) ∈ [0, 2π[×[0, 8π[, with a sufficiently large aspect ratio in order to limit spurious geometrical effects at large shears. The distorted grid is regularly shifted back to the Cartesian grid at times multiple of the shear time 2 3 Ω -1 . Such methods are now standard for local numerical simulations of astrophysical disks [START_REF] Umurhan | Hydrodynamic stability of rotationally supported flows: Linear and nonlinear 2D shearing box results[END_REF][START_REF] Lesur | On the relevance of subcritical hydrodynamic turbulence to accretion disk transport[END_REF] where they are coined as the "shearing box" approach. In our two-dimensional settings, we make use of the vorticity formulation (4), together with the Biot-Savard law to obtain the fluctuating velocity u as a function of the vorticity ω. The stochastic forcing is approximated as shot noise with spatial power spectrum concentrated over wavenumbers k = 4. Time marching uses a second-order Runge-Kutta method, which is explicit for the non-linear term and implicit for the friction and viscous terms.

We have performed several numerical experiments in which we varied the rotation rate Ω while keeping constant the forcing mechanisms, and thus the injection rates ε I and η I . In all cases, transients are followed by the establishment of a statistical steady state associated with a direct cascade of enstrophy. These states are characterized by given values of E, Z, P and R that balance differently in the budgets. Figure 3 shows the terms entering in the energy budget (5) as a function of σ = 3 2 Ωτ f . It provides information about the large scales. The viscous dissipation, which is here proportional to Z, decreases with shear. The contribution from shear through the Reynolds stress instead increases rapidly with σ, and the global balance then requires friction to decrease. The turbulent kinetic energy consequently decreases as a function of the shear rate parameter σ. The terms entering the enstrophy budget (6) are shown in Fig. 4, giving this time information about small-scale turbulent fluctuations. Viscous dissipation increases with σ, because shear stretches vorticity, leading to larger values of the palinstrophy. The enstrophy balance implies that the friction consequently decreases with the rotation rate. The enstrophy Z also decreases and turbulent fluctuations weaken.

Figure 5 reports measurements of the kinetic energy spectrum for the various values of the shear rate parameter σ. The direct cascade of enstrophy in two dimensions is in principle associated to a k -3 scaling range for the spectrum (see, e.g., [START_REF] Boffetta | Two-dimensional turbulence[END_REF]). This law pertains to the case α = 0. However, in the presence of friction, deviations from the -3 scaling have been found in numerics [START_REF] Boffetta | Intermittency in two-dimensional ekman-navier-stokes turbulence[END_REF] and in experiments [START_REF] Boffetta | Effects of friction on 2d turbulence: An experimental study of the direct cascade[END_REF], resulting in a steeper slope of the energy spectrum, with the exponent that decreases as a function of α. Our numerics confirm such findings, with a scaling exponent ≈ -4.6 for σ = 0. As seen from the inset of Fig. 5, the exponents increases with σ and asymptotically reaches -3 for very large shear. This tendency can be explained from previous considerations on the energy and enstrophy budgets (Fig. 3 and Fig. 4) that show a reduction of the strength of friction when σ increases, causing a lower effective value of α. Figure 6 shows the probability density functions (PDF) of the vorticity obtained when varying the shear parameter. The distributions become progressively asymmetrical when increasing σ, with broader negative tails, and reduced positive ones. This is clearly visible from the inset, which shows the skewness ω 3 / ω 2 3/2 becoming more negative and the excess kurtosis ω 4 / ω 2 2 -3 increasing with σ. Shear tends to destroy cyclones (structures with a positive vorticity), favoring the survival of anticyclones. The shear indeed provides a constant negative input of vorticity that entrains anticyclones and damps cyclones. This is at variance with simulations of purely rotating turbulence that favour the persistence of cyclonic vortexes (see, e.g., [START_REF] Biferale | Coherent structures and extreme events in rotating multiphase turbulent flows[END_REF]). Keplerian shear overwhelms this effect. One notices in the figure's inset an unexpected change of tendency at the largest values of σ. We interpret this as geometrical effects that appear at very fast rotations and which could require to consider boxes with an even larger aspect ratio than the one we assumed (4:1).

Above measurements suggest that the role played by anticyclonic vortices gets stronger when shear increases. As seen in the qualitative snapshots of Fig. 2a, b, andc, the shape of such vortices varies with σ. It was found in [START_REF] Kida | Motion of an elliptic vortex in a uniform shear flow[END_REF] that elliptical vortices are possible stable candidates to model such structures. They correspond to an elliptical patch of vorticity, whose major axis is aligned with shear (i.e. with the y direction), and with amplitude

ω = - 3 2 Ω χ + 1 χ (χ -1) , (7)
where χ ≥ 1 denotes the aspect ratio between the major and minor axes. To estimate numerically this aspect ratio, we first identify elliptical points (where the gradient of the total velocity v has two complex-conjugate eigenvalues) corresponding to regions of the flow with a strong-enough rotation. We then compute the vorticity at such locations and invert [START_REF] Hartlep | Cascade model for planetesimal formation by turbulent clustering[END_REF] to express a local value of the aspect ratio χ, which is then averaged in both space and time. Figure 7 shows the dependence of this quantity on the shear rate parameter σ = 3 2 Ωτ f . These results indeed quantify the fact that anticyclones get stretched with increasing shear. Error bars are estimated from partial time averages over 1/5 of the total duration of the simulations. Measurements are very well approximated by a linear fit. When fitting the data, we did not consider the largest value of the rotation rate (σ = 22.8), for which the above procedure leads to a very noisy estimate. Two snapshots of the vorticity are shown to illustrate the typical shape of anticyclones at σ = 0 and σ = 11.4.

III. INERTIAL PARTICLES IN A KEPLERIAN FLOW

We consider tiny dust particles whose dynamics in position-velocity (R p , V p ) phase space follows The three forces appearing in the right-hand side of ( 9) are the linear drag with the gas flow, the Coriolis force, and the combination of gravity and centrifugal forces, respectively.

dR p dt = V p , R p = (x p , y p ) T , (8)
dV p dt = - 1 τ p [V p -v(R p , t)]-2Ω×V p +3Ω 2 x p e x . (9)
The drag is proportional to the difference between the particle velocity and the total gas flow. It involves the particle response time τ p , whose expression depends on the considered regime. If the particle radius a is large-enough compared to the gas mean-free path, their drag follows Stokes regime and τ p = (2/9)ρ p a 2 /(ρ g ν), where ρ p denotes the material mass density of dust grains. Smaller radii follow the Epstein regime where τ p = (8/π) 1/2 ρ p a/(ρ g c s ), with c s denoting the speed of sound. Neglecting both particle-particle interactions and their feedback onto the gas means that we do not need to specify the particle size, making our results relevant to both regimes. The dependence upon the particle response time is usually expressed in terms of a Stokes number obtained by non-dimensionalizing τ p with a characteristic time of the flow. The usual choice in astrophysics is to introduce St Ω = τ p Ω. This Stokes number is mostly relevant at large scales, because it compares the particle response time to the rotation rate, which corresponds to the largest timescale of the flow. Another choice consists in defining St f = τ p /τ f = η 1/3 I τ p . As explained in the previous subsection, the forcing time τ f is actually the only timescale that is fixed in our simulation protocol. Finally, we will mostly rely on the small-scale Stokes number St ω = τ p /τ ω = ω 2 1/2 τ p , which involves the typical amplitude of turbulent fluctuations. Remember that in our protocol, the enstrophy decreases with Ω, so that for a fixed response time, St ω decreases with shear.

As assumed in the previous section, the gas moves in average with the Keplerian speed. Particles which follow (9) do so, and there is hence no mean drag that would lead to a radial drift. Indeed, if we assume that there are no turbulent fluctuations, the gas velocity simply reads v(R p , t) = - 3 2 Ω x p e y and one can easily see that particles moving with the Keplerian shear define stationary solutions to the dynamics. The presence of turbulent fluctuations make particles deviate from the gas streamlines and possibly concentrate. A particular role is played by rotating coherent structures. In usual hydrodynamic situations, vortices act as small centrifuges and normally eject heavy inertial particles. In the present case, the Coriolis force can actually counteract this effect and can sometimes even prevail over the local drag responsible for ejection. This mechanism can be qualitatively understood in terms of Maxey's approximation for the particle velocity [START_REF] Maxey | The gravitational settling of aerosol particles in homogeneous turbulence and random flow fields[END_REF]. At small values of τ p , and in the presence of the Coriolis force, the particle velocity can be approximated as V p ≈ vτ p (dv/dt + 2Ω×v). When Ω = 0, the term ∝ τ p is responsible for driving particles out of vortices. It is however well-known that in the presence of rotation, this term can actually change sign if the Coriolis force overwhelms the fluid acceleration. Particles are then attracted towards the center of the vortices. Such an effect requires that dv/dt and Ω × v have opposite directions, so that it only occurs in anticyclonic vortices.

Previous qualitative argument suggests that particle concentration in anticyclones requires a strong-enough rotation. However, it does not tell anything about dependence on their Stokes number. To get further insight, and building on the results of the previous subsection, we next examine how particles cluster in an elliptical patch of vorticity. We follow again the normal form proposed in [START_REF] Kida | Motion of an elliptic vortex in a uniform shear flow[END_REF] and assume that, locally inside a vortex centered at the origin, the fluid velocity can be approximated as

v = A x y , with A = 3 2 Ω χ -1 0 χ -1 -χ 0 . ( 10
)
As previously noticed, the aspect ratio χ itself depends on the rotation rate Ω. We hereafter use the linear fit proposed in Fig. 7. The next step consists in analyzing the stability of the linear system d dt

R p V p = 0 1 1 τp A +3Ω 2 e xx -1 τp 1 -2ΩJ R p V p , (11)
where e xx = [1 0; 0 0] and J = [0 -1; 1 0]. The vortex center is stable only if all the eigenvalues of the tangent matrix appearing in [START_REF] Umurhan | Streaming instability in turbulent protoplanetary disks[END_REF] have a negative real part.

As shown in Fig. 8, this fixed point is always unstable at the smallest values of Ω (below the black solid line), confirming the idea that vortices eject particles. Above this line, vortex centers are found to always be stable fixed points of the particle dynamics. The convergence rate towards the origin however depends on both the particle response time and the rotation rate. It can be estimated from the real part of the less contracting eigenvalue µ 1 , which is shown as a colored background in Fig. 8. The parameters used in this figure were chosen to match our direct numerical simulations. For instance, we have used the fit of Fig. 4 to express τ ω as a function of Ω. The value of the contraction rate does not give much information, unless it is compared to another characteristic time of the dynamics. Indeed, if we assume that the considered vortex has only a finite lifetime τ c , convergence to the vortex center can only occur if the contraction rate is large enough, namely -Re µ 1 > 1/τ c (illustrated as a red solid line in the figure). Hence, for a fixed value of the shear parameter, only an interval of Stokes numbers can actually lead to strong clustering. At small Stokes numbers, inertia is too weak and particles are migrating too slowly toward the vortex center, namely on timescales that exceed the lifetimes of structures. At large values, particles spiral many times around the origin, making again convergence too slow. As we will now see, such considerations are of relevance to describe clustering in actual turbulent flows.

IV. EJECTION VS. CONCENTRATION IN VORTICES

We now turn to detect and characterize in our simulations the effects of coherent vortices on the particle distribution. As reported in Sec. II, we have performed 18 different sets of simulations where we varied the rotation rate Ω, while keeping constant the forcing characteristics and all other flow parameters. In each simulation, 24 different families of Lagrangian particles with different response times were integrated. They are passive (no feedback onto the flow) and do not interact with each other. rotation (Fig. 9a), particle positions are clearly expelled from vortices, independently of their sign, and concentrate along stripes in between them. At the larger, intermediate values of rotation (σ = 5.7, Fig. 9b), particles are still avoiding vortex cores but distribute this time in spirals around the anticyclones (with a negative vorticity, in blue). Finally, for the largest rotation rate (σ = 14.2, Fig. 9c), all particles have converged in a unique point cluster located in the core of an anticyclone. Increasing the rotation rate has a drastic impact on whether coherent vortices act as centrifuges or give birth to point clusters.

Such observations can be quantified by measuring fluid velocity characteristics along particle paths. The persistent structures of the flow are often linked to its local structure through the Okubo-Weiss parameter, OW (see, e.g., [START_REF] Perlekar | Persistence problem in two-dimensional fluid turbulence[END_REF]). It is defined from the (possibly imaginary) eigenvalues ±β of the strain tensor A ij = ∂ j v i , namely

OW = -β 2 = -(∂ x v x ) 2 -∂ x v y ∂ y v x (12)
This parameter can change sign and therefore divides the flow between OW < 0 (straining, hyperbolic regions, where the two opposite eigenvalues are real) and OW > 0 (rotating, elliptic regions, where the strain has purely imaginary eigenvalues). The Okubo-Weiss parameter has been interpolated at particle positions. Figure 10 shows its average value OW (R p (t), t) over all particles and times, as a function of the shear rate and for various values of the Stokes number. Independently of St ω , this average is negative at small values of σ, meaning that inertial particles spend most of the time in the straining, hyperbolic regions of the flow. However, the average Okubo-Weiss parameter becomes positive for σ 2 and this threshold seems to depend only weakly on the Stokes number. This is consistent with the linear analysis of Kida's vortices from the previous section, which predicts that anticyclones become stable for large-enough shears. At even larger rotation rates, the average Okubo-Weiss parameter shows a fast increase, especially at small values of St ω . Dust particles are therefore living mostly inside rotating, elliptic regions of the flow.

A complementary information on the nature of these eddies can be obtained from the statistics of the fluid vorticity at particle position ω(R p (t), t). The mean and variance over all particles and times, are reported in Fig. 11a andb, as a function of the Stokes number and for various rotation rates. In the absence of rotation (σ = 0), the average ω(R p (t), t) almost vanishes for all Stokes numbers, giving an estimate of the statistical convergence of our data. In the presence of rotation, ω(R p (t), t) is always negative, and this bias increases with σ. Significantly negative average vorticities indicate that the particles spend most of the time in the anticyclones. This happens for fast-enough rotations and intermediate values of St ω . It is for such parameter values that anticyclonic vortices are expected to concentrate particle trajectories, as anticipated in the previous section from the analysis of Kida vortices. The variance of vorticity at particle position is shown in Fig. 11b. At small rotation rates, particles sample flow regions with lesser vorticity fluctu- The PDFs are computed both for the fluid and at the particle position. At low rotation (σ = 1.4, Fig. 12a), both distributions are almost symmetric. Lagrangian statistics display narrower tails than the Eulerian ones, indicating that particle are preferentially sampling the less active regions of the flow. At the larger rotation (σ = 11.4, Fig. 12b), both distributions are skewed toward negative values. However, Lagrangian statistics display a broader negative tail than the Eulerian ones, indicating this time that particles are spending more time in the most active regions of the flow. Skewness toward negative values moreover tells that particles concentrate in anticyclonic vortices, as it could already be seen qualitatively in Fig. 9c.

Despite its apparently noisy aspect, the Lagrangian vorticity distribution shown for St ω = 0.19 and σ = 11.4 in Fig. 12b has actually required a tremendous amount of statistics. For these specific values of the parameters, we have actually performed, in addition to time averages, ensemble average over 20 different realisations of the dynamics. The convergence of all particles in the core of a single, long-living, structure makes average with respect to particles meaningless. For that reason, the results on the mean and variance of the Okubo-Weiss parameter and of the vorticity at particle position (in Fig. 10 and Fig. 11a andb) have been displayed only up to σ ≈ 5. Simulations at larger values of the rotation rate led to much noisier results. As we will see in the next section, the concentration of particles in point clusters can actually be quantified in terms of other observables whose statistics are not affected by such pitfalls.

V. FRACTAL MASS DISTRIBUTION

Results from the previous section indicate that for large-enough rotation rates and intermediate values of their Stokes number, particles preferentially sample the rotating, coherent structures of the flow. This is in line with the linear stability analysis of model vortices that is made in Sec. III. This local analysis, which of course does not straightforwardly apply to fluctuating, time-varying flows, can still be generalized by considering the tangent system to particles' dynamics. The four-dimensional phase-space separation ∆ t = (δR p (t), δV p (t)) T between two infinitesimally close particle trajectories follows

d dt ∆ t = M t ∆ t , (13)
with

M t = 0 1 1 τp A(t) +3Ω 2 e xx -1 τp 1 -2ΩJ , ( 14
) where e xx = [1 0; 0 0], J = [0 -1; 1 0], and A is the gas velocity gradient tensor computed along a reference particle trajectory, namely A ij (t) = ∂ j v i (R p (t), t). The linear evolution 13 integrates to

∆ t = J t ∆ 0 , with J t = T exp t 0 dτ M τ . ( 15
)
J t is the Green operator associated to the linearised dynamics, with T exp the time-ordered exponential. The matrix J T t J t has positive eigenvalues that can be written in the form e 2µj (t)t with j = 1..4. The exponents µ j are called the stretching rates. The Lyapunov exponents are then obtained from the ordered long-time behaviour of the stretching rates: λ j = lim t→∞ µ j (t), with

λ 1 ≥ λ 2 ≥ λ 3 ≥ λ 4 .
The largest Lyapunov exponent λ 1 measures the growth rate of distances between infinitesimally close trajectories. When λ 1 > 0 the dynamics is chaotic, while for λ 1 < 0 phase-space volumes contract and all particle trajectories tend toward each other. Higher-order Lyapunov exponents generalise these notions to infinitesimal phase-space surfaces, volumes, etc. that can either expand or contract. For dissipative dynamics, they can be used to define the Lyapunov dimension [START_REF] Kaplan | Chaotic behavior of multidimensional difference equations[END_REF] as

D L = J - λ 1 + • • • + λ J λ J+1 ( 16
)
where the integer J is such as

λ 1 + • • • + λ J ≥ 0 and λ 1 + • • • + λ J+1 < 0.
This quantity estimates the dimension of the attractor in phase space and has been already successfully used to quantify fractal concentrations of inertial particles [START_REF] Bec | Fractal clustering of inertial particles in random flows[END_REF][START_REF] Boffetta | Large scale inhomogeneity of inertial particles in turbulent flows[END_REF][START_REF] Gustavsson | Ergodic and non-ergodic clustering of inertial particles[END_REF]. Depending on D L , three different regimes for the distribution of particles can be defined. If D L > d, the attractor has a dimension larger than physical space. As particle positions are obtained by projecting one set onto the other, particles fill in that case the whole space. Instead, if D L < d the particles form a fractal set in position-velocity that keeps its fractal structure upon projection, leading to fractal concentrations. Finally, when D L = 0 particles form point clusters in both phase space and positions. Notice that D L cannot take values between 0 and 1. Indeed, if λ 1 > 0, the attractor has to be at least a line. On the contrary, if λ 1 < 0 the dimension of the attractor has to be 0. We have evaluated the Lyapunov exponents in our nu-merical simulations by integrating the tangent dynamics (13) along a large set of particle trajectories. By computing the stretching rates µ j (t) with the standard method of [START_REF] Benettin | Lyapunov characteristic exponents for smooth dynamical systems and for hamiltonian systems; a method for computing all of them. part 1[END_REF], we obtain the Lyapunov exponents as their average over time and particles. Figure 13 represents how the largest Lyapunov exponent λ 1 (expressed in units of τ -1 ω ) depends on both the shear parameter σ and the particle Stokes number St ω . One observes that at large-enough values of the particle response time, τ ω λ 1 decreases with the rotation rate. This depletion of chaoticity mainly comes from the fact that in this limit, turbulent fluctuations get depleted by shear, as already observed in Sec. II. Another noticeable behaviour can be observed at sufficiently large values of σ and small to intermediate St ω (top-left region). We find there that λ 1 < 0 and particle dynamics looses its chaoticity. Figure 15 shows the Lyapunov dimension, for various values of the shear rate parameter σ, this time as a function of the Stokes number. We compare the two different definitions introduced in Sec. III, namely St ω = τ p /τ ω (panel a) and St Ω = τ p Ω (panel b). At slow rotations, small particles form fractal clusters, while larger ones fill space. At larger rotation rates, there is no space-filling regime and the Lyapunov dimension always remains below d = 2. The "turbulent" Stokes number St ω seems the most relevant to describe either small rotation rates or moderate particle response times. For instance, at small to moderate values of σ, the strongest clustering (where D L attains its minimum) is obtained for St ω ≈ 1. Conversely, the asymptotics of large rotation rates and large particle response times seem to depend on St Ω only. Indeed, the bottom-right curves shown in Fig. 15b almost collapse on the top of each other. In the limit σ → ∞, turbulence becomes weaker and the turnover time τ ω increases. The relevant physics is then obtained by balancing particle inertia to the mean shear, whence a dependence upon St Ω only. Such considerations could certainly be formalised in terms of fast and slow timescales allowing one to treat this asymptotics using averaging techniques (for example, through a WKB expansion).

In order to further characterize the particle distribution, we have also studied higher-order statistics of their mass. Standard box-counting techniques rely on the mass (or number) of particles m contained in a box of size . The fractal nature of the particle distribution implies that the moments of m behave as power laws: m n ∼ n , where the exponents n characterize the anomalous scaling typical of multifractal sets. More specifically, n = (n -1)D n + 2 where the D n 's define the dimension spectrum of the attractor. They estimate how strongly a set of n particles clusters.

We first focus on the correlation dimension D 2 . It measures how to the probability p 2 (r) that two particles on the same attractor have a separation distance < depends on scale: Namely we have p 2 (r) ∼ D2 . When D 2 < 2, pairs of particles are more strongly concentrated than those uniformly distributed. We calculate D 2 using box-counting and measure the scaling exponent of m 2 from the average local slope (logarithmic derivative) over Figure 17 shows the fractal dimensions D n , as a function of the order n, for two values of the Stokes number St ω ≈ 0.2 (a) and ≈ 2 (b) and for the same three values of σ. In the case when particles form point clusters (Fig. 17a, σ = 11.4), all fractal dimensions are equal to D n = 0. In the other cases, one clearly observes multifractal clustering with fractal dimensions D n that decrease with n, likely to tend to a limiting value D ∞ at large orders. This behaviour is particularly interesting in the context of the gravitational collapse of dust. Large groups of particles would indeed interact as if distributed on a set of dimension D ∞ , significantly lower that the actual space dimension. Figure 18 shows D 10 , which gives a fair approximation of D ∞ , for the various values of σ. The qualitative shape of variations as a function of St ω is very similar to that of D 2 shown in Fig. 16. Still, despite large error bars, these measurements show that large-order dimensions can be very low over a rather wide range of Stokes numbers. In the absence of rotation (σ = 0), one observes D ∞ ≈ 0.7 at Stokes numbers for which clustering is maximal. Additionally, saturation to D ∞ = 2 seems to occur at Stokes numbers significantly larger than those we considered here. Properties of the mass distribution can be equivalently described in terms of its large deviations. A multifractal distribution implies that the coarse-grained mass m scales as h on a fractal set of dimension D(h), that is with a probability ∝ 2-D(h) . The scaling exponent h is usually known as the singularity exponent, and D(h) is called the dimension spectrum. It is related to the fractal dimensions by a Legendre transform (see, e.g., [START_REF] Paladin | Anomalous scaling laws in multifractal objects[END_REF]): 19 shows the dimension spectrum estimated from the probability distributions of m at various scales . All curves seem defined only for h > h min ≈ 0.7 and D(h) has a finite limit when h → h min . Using the Legendre transform that relates D(h) and D n , one shows that this behavior is equivalent to having the fractal dimension saturating at large order to D ∞ = h min . The values we obtain from these two viewpoints are compatible with each other.

D n+1 = inf h [hn + 2 -D(h)]. Figure

VI. TIME EVOLUTION OF POINT CLUSTERS

We have seen in the previous section that particle distributions develop a non-trivial, multifractal nature that can be used to quantify mass fluctuations. Clearly, this only applies to parameter values that fall in the "fractal clustering" phase of Fig. 14. For values of the Stokes number and of the shear-rate parameter corresponding to the particles concentrating on point clusters, the mass distribution becomes trivial at long times: All fractal dimensions vanish and, equivalently, singularity exponents are restricted to h = 0 with the associated dimension spectrum D(0) = 0. Nevertheless, in this point-clustering regime all the relevant physics occurs during transients.

The Lyapunov exponent λ 1 that we measured in Fig. 13 gives an estimate of how fast a volume of particles converge to the center of a vortex. This tells us about the timescale of local concentration, but the convergence to a unique point cluster is actually ruled by other non-local mechanisms. Our simulations show that, before reaching a steady state (illustrated for instance in Fig. 9c), particles tend to concentrate in several anticyclones. These clusters then survive for long times (of the order of the vortex lifetime) before successively merging together and ending up in the same coherent structure. Such an evolution is illustrated in Fig. 20, which shows how the positions of the various particle clusters (here projected on the x axis) evolve as a function of time. To define clusters, we have here divided the periodic domain [0, 2π[×[0, 8π[ in 256×1024 square boxes. Those containing more than 10% of the total number of particles were then identified as clusters. In the specific case shown in Fig. 20, four different clusters are found to be already formed after t ≈ 40τ ω . We then track their time evolution and observe that they successively merge and finally form a unique point cluster at the end of the simulation. The typical rate at which clusters merge together is given by the lifetime of vortices. This time is found to be τ c ≈ 225 τ ω . It is much longer than the time -1/λ 1 ≈ 2.3 τ ω associated to particle concentration inside a given vortex core, but it is also much longer than all the physically relevant timescales of our problem. This separation of timescales suggests that cluster merger can actually be well approximated as a sequence of Poisson processes. Let us denote by N (t) the average number of clusters that are present in the domain at time t, where we assume that particles are uniformly seeded in the flow at t = 0. Above considerations suggest that once clusters are formed, their average number approximately follows the kinetics

d N dt ≈ -γ N ( N -1), (17)
where the rate γ is of the order of τ -1 c . The above dynamics just assumes that between times t and t + dt, all pairs of clusters can independently merge with a probability γ dt. Equation [START_REF] Tanga | Forming planetesimals in vortices[END_REF]

integrates to N ≈ 1 1 -c e -γ t . (18)
We have repeated the cluster-tracking procedure that is explained above over an ensemble of 24 different realisations with the same physical parameters (σ = 11.4 and St ω = 0.19). Figure 21 shows the time evolution of the average number of clusters then obtained.

Data are well fitted by formula [START_REF] Bracco | Particle aggregation in a turbulent keplerian flow[END_REF] with c = 0.81 and λ = 0.0014 τ -1 ω ≈ 0.3 τ -1 c (dashed line). As anticipated, this rate seems ruled by the lifetime of anticyclones. However, because such measurements require a large number of realisations, we have not investigated how these numbers depend upon the shear parameter σ. Nevertheless, we expect stronger shears to deplete the survival of coherent structures, as for instance reported in [START_REF] Campana | Stochastic model for the alignment and tumbling of rigid fibres in two-dimensional turbulent shear flow[END_REF]. This implies that an increasing rotation rate could make the merger of point clusters a faster process.

VII. CONCLUDING REMARKS

In this paper we studied the clustering of heavy, dust particles in a turbulent flow with Keplerian rotation and shear. To do so, we performed 2D direct numerical simulations. The results have been analysed with tools borrowed from the study of nonlinear dynamical systems. The main findings are:

-The turbulent flow is substantially modified when increasing shear. Strong rotation favours anisotropies, with the preferential formation and survival of anticyclonic vortices.

-Particles form point clusters in anticyclones or concentrate on fractal sets depending on both rotation and particle response time.

-The mass distribution of dust particles is multifractal, with dimensions that saturate at large orders, and decrease with rotation.

-Clusters form in a hierarchical manner and merge one with the other in the course of time. This evo-lution occurs on timescales much longer than all other physical processes.

These results show that turbulence in protoplanetary disks is able to concentrate particles of various sizes. In the presence of fast rotation, these particles are rapidly and strongly clustered inside anticyclones. We expect these clusters to easily form planetesimals by gravitational collapse. We believe that planetesimals could form also from fractal sets, due to the important clustering strength between large groups of particles. We give below more details on the importance and limits of our results on the astrophysical problem that was our primary motivation.

A first remark can be drawn from our extensive use of tools borrowed from nonlinear physics. This represents a rather powerful approach to tackle problems related to astrophysical disks dynamics. The new analysis tools that we could identify give original, innovative viewpoints on planetesimal formation.

The effectiveness of such analysis tools does not depend on the details of the physical environment that we have considered. They remain of relevance in more complex situations. We indeed expect many aspects of dust clustering (point clustering in anticyclones, multifractal distributions, saturation of high-order dimensions, etc.) to be universal and independent of the turbulencegenerating mechanisms (vertical shear instability, MRI, . . . ). In our model we considered a generic way to maintain a turbulent state, but different forcing terms could have been used to mimic the effects of specific instabilities. This would make our results more easily comparable with previous studies in astrophysics and using our tools to analyse such results could clearly provide new insights.

In our approach we used a fixed, white-noise forcing, thereby prescribing the inputs of energy and enstrophy to our turbulent flow. This allowed us to identify the key physical parameters for dust concentration, namely the rotation rate of the flow Ω and the particle response time to the flow dynamics τ p . Most studies of planetesimal formation consider a fixed rotation rate. At variance we have here changed it in a systematic manner in order to highlight how the combination of rotation and shear affect the properties of both the flow and the concentration of solid particles.

Performing such a systematic study of course required to simplify the physical setup and to neglect several mechanisms at play in protoplanetary disks. First of all, we considered the gas to be in Keplerian rotation, as the dust particles are. With such a simplification, solids do not drift radially toward the star. Actually this average motion would have added a source of gas/dust velocity difference and thus modified the clustering of particles in anticyclones. Nevertheless, the drift timescale is known to be orders of magnitude larger than the rotation one, allowing to decouple the two processes. Another simplification was to neglect the back reaction of the dust onto the gas. Again, this is a second-order effect as the particles are strongly diluted. However, considering this coupling is needed to describe other planetesimal-forming processes, in particular the streaming instability which could be triggered locally in highly-concentrated regions. Finally, self-interactions between dust particles (e.g. collisions, gravity, etc.) have been disregarded. They should of course be included to follow dynamical evolution of particles after clustering. All these various aspects are kept for future work.

As discussed in subsection II A, we focused in this study on two-dimensional flows. Protoplanetary disks, despite being thin, have a vertical extension that is of great relevance for dust stratification. Solid particles indeed settle toward the mid-plane, on a timescale varying with the particle size. Numerous instabilities can also develop along the vertical direction, resulting in a turbulent state that can stir up again some solids to the disk surface. Full three-dimensional simulations should therefore be performed to address the question of the effect of vertical stratification on clustering. We nevertheless highlight again that at high rotation rate (Ω 1) the flow is mainly two-dimensional, due to the formation of Taylor columns. We thus expect the formation of strong clusters of particles (lines in the core of columnar vortices) to occur also in three dimensions, at least for large-enough rotation frequencies. Moreover, new disk observations [START_REF] Villenave | A highly settled disk around oph163131[END_REF] suggest that the vertical extent of dust is even smaller than previously thought. This requires solids to settle almost completely to the mid-plane, making their evolution nearly 2D from that moment on.

Let us finally stress that protoplanetary disks are compressible, with the propagation of spiral density waves at large scales. However, as already mentioned, we consider here sufficiently small scales to justify the incompressible approximation (see for instance [START_REF] Meheut | Angular momentum transport and large eddy simulations in magnetorotational turbulence: the small Pm limit[END_REF] for a comparison between compressible and incompressible turbulence in the shearing box). Addressing larger scales will of course require accounting for compressibility.

FIG. 2 .

 2 FIG. 2. Three different snapshots of the gas vorticity field obtained in the absence of rotation and shear (a), for an intermediate value of Ω (b), and for a large value of Ω (c).

FIG. 3 .FIG. 4 .

 34 FIG.3. Contributions of friction, dissipation, and shear to the energy balance (5) as a function of the shear rate parameter σ.

FIG. 5 .

 5 FIG.5. Time-averaged power spectra of the turbulent kinetic energy for various values of the shear parameter σ, as labelled. Inset: Measured scaling exponent that shows a continuous evolution from a value lesser than -4.5 to -3 when the shear rate increases.

FIG. 6 .

 6 FIG. 6. Probability density functions (PDF) of the vorticity for various values of the shear parameter σ. The inset shows the vorticity skewness and excess kurtosis as a function of the shear parameter.

FIG. 7 .

 7 FIG. 7. Aspect ratio χ of anticyclonic vortices as a function of the shear rate parameter σ = 3 2 Ωτ f . The solid line is a linear law fitting the data of the form χ = 1 + 0.19 σ. The two insets show snapshots of the vorticity perturbation ω (normalized by τ -1 ω ) for σ = 0 and σ = 11.4.

FIG. 8 .

 8 FIG.8. Linear analysis of particles dynamics in a stationary vortex patch given by (11): Contour plots of the eigenvalue µ1 with largest real part in the Stokes number / shear rate plane. The black solid line corresponds to Re µ1 = 0, while the red one delimit a region where -Re µ1 is larger than a given threshold (see text).

FIG. 9 .

 9 FIG. 9. Vorticity ω, by τ -1 ω , and dust particle position Rp (black dots) for Stω ≈ 0.35 and (a) σ=0, (b) σ=5.7, and (c) σ=14.2.

2 FIG. 10 .

 210 FIG. 10. Average Okubo-Weiss parameter at particle position OW (Rp) as a function of the shear rate parameter σ for various values of the Stokes number Stω.

7 FIG. 11 .

 711 FIG. 11. Mean (a) and variance (b) of vorticity ω(Rp(t), t) at particle position, as a function of the Stokes number Stω, for various values of the shear rate parameter σ = 3 2 Ωτ f .

FIG. 12 .

 12 FIG. 12. PDFs of the Eulerian vorticity ω and of the Lagrangian one ω(Rp(t), t) for Stω ≈ 0.19 and (a) σ = 1.4, (b) σ = 11.4.

1 FIG. 13 .

 113 FIG. 13. Phase diagrams of the shear rate parameter σ vs the Stokes number Stω. The color code is for the largest Lyapunov exponent λ1 (normalized by τ -1 ω ). Each black dot stands for the parameters that are actually simulated.

FIG. 14 .

 14 FIG. 14. Phase diagrams of the shear rate parameter σ vs the Stokes number Stω. The color code is for the Lyapunov dimension DL. Each black dot stands for the parameters that are actually simulated.

Figure 14

 14 Figure 14 shows the Lyapunov dimension D L as a function of σ and St ω . At low rotation rates and large Stokes numbers, the particles are space-filling and D L > 2. They otherwise form fractal clusters where 0 < D L < 2. At parameters values where the largest Lyapunov exponent is negative, we have D L = 0 and particles concentrate in point clusters. It is clear from this phase diagram that rotation favours the clustering of heavy particles. When comparing it with the predictions from linear analysis of Fig. 8, one recognizes similar behaviours, and in particular point clusters are found to occur at similar parameter values (fast-enough rotation and intermediate Stokes numbers).Figure15shows the Lyapunov dimension, for various values of the shear rate parameter σ, this time as a function of the Stokes number. We compare the two different definitions introduced in Sec. III, namely St ω = τ p /τ ω (panel a) and St Ω = τ p Ω (panel b). At slow rotations, small particles form fractal clusters, while larger ones fill space. At larger rotation rates, there is no space-filling

8 FIG. 15 .

 815 FIG. 15. Lyapunov dimension DL as a function of (a) the "turbulent" Stokes number Stω and (b) the "rotational" Stokes number StΩ, for various values of the shear rate parameter σ as labelled.

4 FIG. 16 . 3 2

 4163 FIG. 16. Correlation dimension D2 as a function of the Stokes number Stω for three values of the shear rate parameter σ = 3 2 Ωτ f .

4 FIG. 17 .

 417 FIG. 17. Fractal dimension Dn as a function of the order of the fractal dimension n for three values of the shear rate parameter σ and Stokes number Stω ≈ 0.2 (a) and Stω ≈ 2 (b).

4 FIG. 18 .

 418 FIG. 18. D∞ as a function of the Stokes number Stω for different values of the shear rate parameter σ.

FIG. 19 .

 19 FIG. 19. Dimension spectrum D(h) as a function of the singularity exponent h for different values of the box size for σ = 3.6 and Stω = 0.16.

4 FIG. 20 .

 420 FIG. 20. Position in the x-direction of particle clusters as a function of time for σ = 11.4 and Stω = 0.19.

FIG. 21 .

 21 FIG. 21. Average number of clusters as a function of time for σ = 11.4 and Stω = 0.19. The dashed line is a fit of the data obtained from the kinetic approximation (18) with λτω = 0.0014 and c = 0.81.

ACKNOWLEDGEMENTS

This work received support from the UCA-JEDI Future Investments at the Université Côte d'Azur (funded by the French government, and managed by the National Research Agency, ANR-15-IDEX-01). This work was supported by the "Programme National de Planétologie" (PNP), "Programme National de Physique Stellaire" (PNPS) and "Programme National de Physique et Chimie du Milieu Interstellaire" (PCMI) of CNRS/INSU co-funded by CEA and CNES. HM acknowledges inspiring discussions that took place at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611. FG and JB would like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality (EPSRC grant No EP/R014604/1) during the programme "Mathematical aspects of turbulence: where do we stand?" where a part of this work was undertaken.


Jaya Kumar
email: jayaka@iisc.ac.in

Akhilesh Kumar Verma
email: akhilesh@iisc.ac.in

Jérémie Bec
email: jeremie.bec@mines-paristech.fr

Rahul Pandit

Machine learning strategies for path-planning microswimmers in turbulent flows

We develop an adversarial-reinforcement learning scheme for microswimmers in statistically homogeneous and isotropic turbulent fluid flows, in both two (2D) and three dimensions (3D). We show that this scheme allows microswimmers to find non-trivial paths, which enable them to reach a target on average in less time than a naïve microswimmer, which tries, at any instant of time and at a given position in space, to swim in the direction of the target. We use pseudospectral direct numerical simulations (DNSs) of the 2D and 3D (incompressible) Navier-Stokes equations to obtain the turbulent flows. We then introduce passive microswimmers that try to swim along a given direction in these flows; the microswimmers do not affect the flow, but they are advected by it. Two, non-dimensional, control parameters play important roles in our learning scheme: (a) the ratio Ṽs of the microswimmer's bare velocity Vs and the root-mean-square (rms) velocity urms of the turbulent fluid; and (b) the product B of the microswimmer-response time B and the rms vorticity ωrms of the fluid. We show that the average time required for the microswimmers to reach the target, by using our adversarial-learning scheme, eventually reduces below the average time taken by microswimmers that follow the naïve strategy.

I. INTRODUCTION

Machine-learning techniques and advances in computational facilities have led to significant improvements in obtaining solutions to optimization problems, e.g., to problems in path planning and optimal transport, referred to in control systems as Zermelo's navigation problem [START_REF] Pandit | Zermelo Über das navigationsproblem bei ruhender oder veränderlicher windverteilung[END_REF]. With vast amounts of data available from experiments and simulations in fluid dynamics, machinelearning techniques are being used to extract information that is useful to control and optimize flows [START_REF] Brunton | Machine Learning for Fluid Mechanics[END_REF]. Recent studies include the use of reinforcement learning, in fluidflow settings, e.g., (a) to optimise the soaring of a glider in thermal currents [START_REF] Reddy | Learning to soar in turbulent environments[END_REF] and (b) the development of an optimal scheme in two-(2D) and three-dimensional (3D) fluid flows that are time independent [START_REF] Colabrese | Flow navigation by smart microswimmers via reinforcement learning[END_REF][START_REF] Gustavsson | Finding efficient swimming strategies in a three-dimensional chaotic flow by reinforcement learning[END_REF]. Optimal locomotion, in response to stimuli, is also important in biological systems ranging from cells and micro-organisms [START_REF] Dusenbery | Living at Micro Scale: The Unexpected Physics of Being Small[END_REF][START_REF] Durham | Turbulence drives microscale patches of motile phytoplankton[END_REF][START_REF] Michalec | Turbulence triggers vigorous swimming but hinders motion strategy in planktonic copepods[END_REF] to birds, animals, and fish [START_REF] Verma | Efficient collective swimming by harnessing vortices through deep reinforcement learning[END_REF]; such locomotion is often termed taxis [10].

It behooves us, therefore, to explore machine-learning strategies for optimal path planning by microswimmers in turbulent fluid flows. We initiate such a study for microswimmers in 2D and 3D turbulent flows. In particular, we consider a dynamic-path-planning problem that seeks to minimize the average time taken by microswimmers to reach a given target, while moving in a turbulent fluid flow that is statistically homogeneous and isotropic. We develop a novel, multi-swimmer, adversarial-Q-learning algorithm to optimise the motion of such microswimmers that try to swim towards a specified target (or tar-gets). Our adversarial-Q-learning approach ensures that the microswimmers perform at least as well as those that adopt the following naïve strategy: at any instant of time and at a given position in space, a naïve microswimmer tries to point in the direction of the target. We examine the efficacy of this approach as a function of the following two dimensionless control parameters: (a) Ṽs = V s /u rms , where the microswimmer's bare velocity is V s and the the turbulent fluid has the root-mean-square velocity u rms ; and (b) B = B ω rms , where B is the microswimmerresponse time and ω rms the rms vorticity of the fluid. We show, by extensive direct numerical simulations (DNSs), that the average time T , required by a microswimmer to reach a target at a fixed distance, is lower, if it uses our adversarial-Q-learning scheme, than if it uses the naïve strategy.

II. BACKGROUND FLOW AND MICROSWIMMER DYNAMICS

For the low-Mach-number flows we consider, the fluidflow velocity u satisfies the incompressible Navier-Stokes (NS) equation. In two dimensions (2D), we write the NS equations in the conventional vorticity-stream-function form, which accounts for incompressibility in 2D [START_REF] Pandit | An overview of the statistical properties of two-dimensional turbulence in fluids with particles, conducting fluids, fluids with polymer additives, binary-fluid mixtures, and superfluids[END_REF]:

(∂ t + u • ∇) ω = ν∇ 2 ω -α ω + F ω ; (1)
here, u ≡ (u x , u y ) is the fluid velocity, ν is the kinematic viscosity, α is the coefficient of friction (present in 2D, e.g., because of air drag or bottom friction) and the vorticity ω = (∇ × u), which is normal to u in 2D. The 3D incompressible NS equations are p is the pressure and the density ρ of the incompressible fluid is taken to be 1; the large-scale forcing F ω (largescale random forcing in 2D) or f (constant energy injection in 3D) maintains the statistically steady, homogeneous, and isotropic turbulence, for which it is natural to use periodic boundary conditions. We consider a collection of N p passive, non-interacting microswimmers in the turbulent flow; X i and pi are the position and swimming direction of the microswimmer. Each microswimmer is assigned a target located at X T i . We are interested in minimizing the time T required by a microswimmer, which is released at a distance r 0 = |X i (0) -X T i | from its target, to approach within a small distance r = |X i (T ) -X T i | r 0 of this target. The microswimmer's position and swimming direction evolve as follows [START_REF] Pedley | Hydrodynamic phenomena in suspensions of swimming microorganisms[END_REF]:

(∂ t + u • ∇) u = -∇p/ρ + f + ν∇ 2 u; ∇.u = 0; (2
dX i dt = u(X i , t) + V s pi ; ( 3
)
dp i dt = 1 2B [ô i -(ô i .p i ) pi ] + 1 2 ω × pi ; (4)
here, we use bi-linear (tri-linear) interpolation in 2D (3D) to determine the fluid velocity u at the microswimmer's position X i from eq. 2; V s pi is the swimming velocity, B is the time-scale associated with the microswimmer to align with the flow, and ôi is the control direction. Equation 4 implies that pi tries to align along ôi . We define the following non-dimensional control parameters: Ṽs = V s /u rms , where u rms = |u| 2 1/2 is the root-meansquare (rms) fluid flow velocity, and B = B/τ Ω , where τ Ω = ω -1 rms ; ω rms = |ω| 2 1/2 denotes the root-meansquare vorticity.

III. ADVERSARIAL Q-LEARNING FOR SMART MICROSWIMMERS

Designing a strategy consists in choosing appropriately the control direction ôi , as a function of the instantaneous state of the microswimmer, in order to minimize the mean arrival time T . To develop a tractable framework for Q-learning, we use a finite number of states by discretizing the fluid vorticity ω at the microswimmer's location into 3 ranges of values labelled by S ω and the angle θ i , between pi and Ti , into 4 ranges S θ , as shown in fig 1. The choice of ôi is then reduced to a map from (S ω , S θ ) to an action set, A, which we also discretize into the following four possible actions: A := Ti , -Ti , Ti⊥ , -Ti⊥ , where Ti = (X T i -X i )/|X T i -X i | is the unit vector pointing from the swimmer to its target and ( Ti⊥ • Ti ) = 0. Therefore, for the naïve strategy ôi (s i ) ≡ Ti , ∀ s i ∈ (S ω , S θ ). This strategy is optimal if Ṽs 1: Microswimmers have an almost ballistic dynamics and move swiftly to the target. For Ṽs 1, vortices affect the microswimmers substantially, so we have to develop a nontrivial Q-learning strategy, in which ôi is a function of ω(X i , t) and θ i . In our Q-learning scheme, we assign a quality value to each state-action binary relation of microswimmer i as follows:

Q i : (s i , a i ) → R, where s i ∈ (S ω , S θ )
and a i ∈ A; and we use the -greedy method [START_REF] Watkins | Learning from Delayed Rewards[END_REF] (with parameter g ), in which the control direction is chosen from the probability distribution P [ô i (s i )] = g /4 + (1g ) δ (ô i (s i ) -ômax ), where ômax := argmax a∈A Q i (s i , a) and δ(.) is the Dirac delta function. At each iteration, ôi is calculated as above and the microswimmer evolution is performed by using eqs. 3 and 4. In the canonical Q-learning approach, during the learning process, each of the Q i 's are evolved by using the Bellman equation [START_REF] Sutton | Reinforcement learning: An introduction[END_REF] below, whenever there is a state change, i.e., s i (t) = s i (t + δt): [START_REF] Gustavsson | Finding efficient swimming strategies in a three-dimensional chaotic flow by reinforcement learning[END_REF] where λ and γ are learning parameters that are set to optimal values after some numerical exploration (see tab. I), and R i is the reward function. For the path-planning problem we define

Q i (s i (t), ôi (s i (t))) → (1 -λ) Q i (s i (t), ôi (s i (t))) + λ R i (t) + γ max a∈A Q i (s i (t + δt), a) ,
R i (t) = |X i (t -n δt) -X T i | -|X i (t) - X T i |,
where n = min l∈N {s i (t -l δt) = s i (t)}. According to eq. 5, any ôi for which R i is positive can be a solution, and there exist many such solutions that are sub-optimal compared to the naïve strategy.

To reduce the solution space, we propose an adversarial scheme: Each microswimmer, the master, is accompanied by a slave microswimmer, with position X Sl i (t), that shares the same target at X T i , and follows the naïve strat-γ = 0.99 λ = 0.01 g = 0.001 ω0/ωrms = 1.0 TABLE I. List of learning parameter values: γ is the earning discount, λ is the learning rate, g is the -greedy algorithm parameter that represents the probability with which the nonoptimal action is chosen, ω0 is the cut-off used for defining Sω, and ωrms is the rms value of ω.

egy, i.e., ôSl

i (t) ≡ TSl i = (X Sl i -X T i )/|X Sl i -X T i |.
Now, whenever the master undergoes a state change, the corresponding slave's position and direction are re-initialized to that of the master, i.e., if s i (t) = s i (t + δt), then X Sl i (t + δt) = X i (t + δt) and pSl i (t + δt) = pi (t + δt) (see fig. 2). Then the reward function for the master microswimmer is given by

R AD i (t) = |X Sl i (t) -X T i | - |X i (t) -X T i |; i.e.
, only those changes that improve on the naïve startegy are favored.

In the conventional Q-learning approach [START_REF] Kaelbling | Reinforcement learning: A survey[END_REF][START_REF] Watkins | Learning from Delayed Rewards[END_REF], the matrices Q i of each microswimmer evolve independently; this matrix is updated only after a state change, so a large number of iterations are required for the convergence of Q i . To speed-up this learning process, we use the following multi-swimmer, parallel-learning scheme: all the microswimmers share a common Q matrix, i.e., Q i = Q, ∀i. At each iteration, we choose one microswimmer at random, from the set of microswimmers that have undergone a state change, to update the corresponding element of the Q matrix (flow chart in Appendix A); this ensures that the Q matrix is updated at almost every iteration and so it converges rapidly.

IV. NUMERICAL SIMULATION

We use a pseudospectral DNS [START_REF] Canuto | Spectral methods[END_REF][START_REF] Pandit | Statistical properties of turbulence: an overview[END_REF], with the 2/3 dealiasing rule to solve eqs. 1 and 2. For time marching we use a third-order Runge-Kutta scheme in 2D and the exponential Adams-Bashforth time-integration scheme in 3D; the time step δt is chosen such that the Courant-Friedrichs-Lewy (CFL) condition is satisfied. Table II gives the parameters for our DNSs in 2D and 3D, such as the number N of collocation points and the Taylormicroscale Reynolds numbers R λ = u rms λ/ν, where the Taylor microscale

λ = k k 2 E(k)/ k E(k) -1/2 .

A. Naïve microswimmers

The average time taken by the microswimmers to reach their targets is T (see fig. 3). If Ti = (X i -X T i )/|X i -X T i | is the unit vector pointing from the microswimmer to the target, then for Ṽs 1 we expect the naïve strategy, i.e., ôi = Ti , to be the optimal one. For Ṽs 1, we observe that the naïve strategy leads (a) Illustrative (blue) paths for two microswimmers, with their corresponding (yellow) circular target regions (mapping in red dashed lines) where the microswimmer is eventually absorbed and re-initialized. We consider random positions of targets and initialize a microswimmer at a fixed distance from its corresponding target with randomized p; (b) a snapshot of the microswimmer distribution, in a vorticity field (ω), for the naïve strategy, at time t = 30τΩ, with Ṽs = 1. Here, the initial distance of the microswimmers from their respective targets is L/3 and the target radius is L/50; we use a system size L with periodic boundary conditions in all directions.

B. Smart microswimmers

In our approach, the random initial positions of the microswimmers ensures that they explore different states without reinitialization for each epoch. Hence, we present results with 10000 microswimmers, for a single epoch. In our single-epoch approach, the control map ôi reaches a steady state once the learning process is complete (fig. 5(b)). We would like to clarify here that, in our study, the training is performed in the fully turbulent time-dependent flow; even though this is more difficult than training in a temporally frozen flow, the gains, relative to the naïve strategy, justify this additional level of difficulty.

We use the adversarial Q-learning approach outlined above (parameter values in tab. I) to arrive at the optimal scheme for path-planning in a 2D turbulent flow. To quantify the performance of the smart microswimmers, we introduce equal numbers of smart (master-slave pairs) and naïve microswimmers into the flow. The scheme presented here pits Q-learning against the naïve strategy and enables the adversarial algorithm to find a strategy that can out-perform the naïve one. (Without the adversarial approach, the final strategy that is obtained may end up being sub-optimal.) FIG. 5.

Learning statistics: (a) Plot of T |t, ∆ , with ∆ = 10 τΩ, in 2D. Adversarial Q-learning initially shows a transient behavior, before settling to a lower value of T than that in the naïve strategy. (b) The evolution of the control map, ôi, where the color codes represent the actions that are performed for each of the 12 states. Initially, Qlearning explores different strategies and settles down to a ôi that shows, consistently, improved performance relative to the naïve strategy. FIG. 6. The dependence of T on Ṽs, for different values of B, shown for the naïve strategy (dotted line) and for adversarial Q-learning (solid line), for our 2D turbulent flow. The plot shows that, in the parameter space that we have explored, our adversarial-Q-learning method yields a lower value T than in the naïve strategy. The plot in the inset shows that the CPDF of T has an exponential tail.

V. RESULTS

The elements of Q evolve during the initial-training stage, so P > (T ) also evolves in time until the system reaches a statistically steady state (in which the elements of Q do not change). Hence, T also changes during the initial-training stage; to capture this time dependence, we define T (t) := 1/N (t)

N (t)
i=1 T i , where T i is the time taken by the i th microswimmer, since its intialization, to arrive at its target at the time instant t and N (t) is the number of microswimmers that reach their targets at time instant t. We find that T (t) shows large fluctuations; so we average it over a time window ∆ and define T |t, ∆ : Figures 5(a), and 5(b) show the evolution of T |t, ∆ and ô, respectively, for the naïve strategy and our adversarial-Q-learning scheme. After the initial learning phase, the Q-learning algorithm explores different ô, before it settles down to a steady state. It is not obvious, a priori, if there exists a stable, non-trivial, optimal strategy, for microswimmers in turbulent flows, that could out-perform the naïve strategy. The plot in fig. 6 shows the improved performance of our adversarial-Q-learning scheme over the naïve strategy, for different values of Ṽs and B; in these plots we use T = T |t → ∞, ∆ , so that the initial transient behavior in learning is excluded. The inset in fig. 6 shows that P > (T ) has an exponential tail, just like the naïve scheme in fig. 4, which implies the smart microswimmers also get trapped; but a lower value of T implies they are able to escape from the traps faster than microswimmers that employ the naïve strategy. Note that the presence of a possible noise in the measurement of the discrete vorticity S ω should not change our findings because of the coarse discretization we use in defining the states.

In a 3D turbulent flow, we also obtain such an improvement, with our adversarial Q-learning approach, over the naïve strategy. The details about the 3D flows, parameters, and the definitions of states and actions are given in Appendix B. In fig. 7 we show a representative plot, for the performance measure, which demonstrates this improvement in the 3D case (cf. fig. 5 for a 2D turbulent flow). The evolution of ôi in 3D shows that learning has not reached a steady state due to lower probability of swimmers reaching the target, compared to 2D case.

VI. CONCLUSIONS

We have shown that the generic Q-learning approach can be adopted to solve control problems arising in com-plex dynamical systems. In [START_REF] Kularatne | Optimal path planning in time-varying flows using adaptive discretization[END_REF], global information of the flows has been used for path-planning problems in autonomous-underwater-vehicles navigation to improve their efficiency, based on the Hamilton-Jacobi-Bellmann approach. In contrast, we present a scheme that uses only the local flow parameters for the path planning.

The flow parameters (tab. II) and the learning parameters (tab. I) have a significant impact on the performance of our adversarial-Q-learning method. Even the choice of observables that we use to define the states (S ω , S θ ) can be changed and experimented with. Furthermore, the discretization process can be eliminated by using deeplearning approaches, which can handle continuous inputs and outputs [START_REF] Lillicrap | Continuous control with deep reinforcement learning[END_REF]. Our formulation of the optimal-pathplanning problem for microswimmers in a turbulent flow is a natural starting point for detailed studies of control problems in turbulent flows.

VII. DISCUSSION

We were made aware of [START_REF] Biferale | Zermelo's problem: Optimal point-topoint navigation in 2D turbulent flows using reinforcement learning[END_REF] during the writing of this manuscript, where they tackle the problem using an Actor-Critic reinforcement learning scheme.

We contrast, below, our reinforcement-learning approach with that of Ref. [START_REF] Biferale | Zermelo's problem: Optimal point-topoint navigation in 2D turbulent flows using reinforcement learning[END_REF].

• Reference [START_REF] Biferale | Zermelo's problem: Optimal point-topoint navigation in 2D turbulent flows using reinforcement learning[END_REF] uses 900 discrete states, which are defined based on the approximate location of the microswimmer. By contrast, our scheme uses only the local vorticity (S ω ), at the position of the microswimmer, and the orientation (S θ ); after discretization, we retain only 12 states. In analogy with navigation parlance, Ref. [START_REF] Biferale | Zermelo's problem: Optimal point-topoint navigation in 2D turbulent flows using reinforcement learning[END_REF] uses a GPS and our approach uses a light-house along with a localvorticity measurement.

• In Ref. [START_REF] Biferale | Zermelo's problem: Optimal point-topoint navigation in 2D turbulent flows using reinforcement learning[END_REF], the states are sensed periodically and the elements of Q are updated at every sensing instant. In constrast, we monitor the states continuously and update the elements of Q only when there is a state change. If the periodicity of sensing is smaller than the rate of change in states of the microswimmer, both schemes should show similar convergence behaviors.

• Reference [START_REF] Biferale | Zermelo's problem: Optimal point-topoint navigation in 2D turbulent flows using reinforcement learning[END_REF] uses a conventional, episode-based training approach, which is sequential, whereas we use multiple microswimmers to perform parallel training.

• Reference [START_REF] Biferale | Zermelo's problem: Optimal point-topoint navigation in 2D turbulent flows using reinforcement learning[END_REF] uses an actor-critic approach, whereas we use an adversarial-learning method.

  ) arXiv:1910.01728v2 [physics.flu-dyn] 7 May 2021

FIG. 1 .

 1 FIG. 1. Left panel: a pseudocolor plot of the vorticity field, with a microswimmer represented by a small white circle; the black arrow on the microswimmer indicates its swimming direction, p, the red arrow represents the direction towards the target, T, and θ is the angle between p and T. Top-center panel shows the discretized vorticity states (red |||: ω > ω0, green \\\: -ω0 ≤ ω ≤ ω0, blue ///: ω < -ω0). In our approach we use ω0 = ωrms. The bottom-center panel indicates the color code for the discretized θ(red |||: -π/4 ≤ θ < π/4; orange \\\: π/4 ≤ θ < 3π/4; blue ///: -3π/4 ≤ θ < -π/4; gray ≡: 3π/4 ≤ θ < 5π/4). The right panel lists all possible discrete states of the microswimmers, via colored squares where the lower half stands for the vorticity state, Sω, and the upper half represents the direction state, S θ .

FIG. 2 .

 2 FIG. 2. Top-left panel: a schematic diagram illustrating the trajectories of master (black line) and slave (dashed black line) microswimmers superimposed on a pseudocolor plot of the two-dimensional (2D) discrete vorticity field Sω; the master undergoes a state change at the points shown by white filled circles; white arrows indicate the re-setting of the slave's trajectory. Top-right panel: color code for the control direction ôi; for the states si ∈ (Sω, S θ ) see Fig. 1. Bottom panel: control maps for the master and slave; for the purpose of illustration, we use ôi = Ti⊥ , for the master; for Ṽs 1 and B = 0, this leads to the circular path shown in our schematic diagram.

FIG. 3 .

 3 FIG. 3.(a) Illustrative (blue) paths for two microswimmers, with their corresponding (yellow) circular target regions (mapping in red dashed lines) where the microswimmer is eventually absorbed and re-initialized. We consider random positions of targets and initialize a microswimmer at a fixed distance from its corresponding target with randomized p; (b) a snapshot of the microswimmer distribution, in a vorticity field (ω), for the naïve strategy, at time t = 30τΩ, with Ṽs = 1. Here, the initial distance of the microswimmers from their respective targets is L/3 and the target radius is L/50; we use a system size L with periodic boundary conditions in all directions.

FIG. 4 .

 4 FIG. 4. Plots showing exponential tails in P > (T ) for the naïve strategy, with different values of Ṽs and B. The inset shows how these data collapse when, T is normalized, for each curve, by the corresponding T , which implies P > (T ) ∼ exp (-T / T ).

  = 1/∆ t+∆ t T (τ ) dτ . The initial growth in T |t, ∆ arises because T |t, ∆ ≤ t. The plots in figs. 5(a) and 7 show the time evolution of T |t, ∆ for the smart and naïve microswimmers. Note that ôi becomes a constant, for large t, in fig. 5(b); this implies that the elements of Q have settled down to their steady-state values.

FIG. 7 .

 7 FIG. 7. Learning statistics in 3D: (a) The performance trend, T |t, ∆ /τΩ, with ∆ = 10τΩ for adversarial Q-learning (blue line) and naïve strategy (red broken-line) for microswimmers in a 3D homogeneous isotropic turbulent flow, for Ṽs = 1.5 and B = 0.5. The trend shows a slow rise in performance, similar to that observed in 2D. In 3D the Q-learning is performed by using 13 states and 6 actions defined in Appendix B; (b) The evolution of ôi in 3D shows that learning has not reached a steady state due to lower probability of swimmers reaching the target, compared to 2D case.

ACKNOWLEDGMENTS

We thank DST, CSIR (India), BRNS, and the Indo-French Centre for Applied Mathematics (IFCAM) for support.

Appendix A: Flowchart Figure 8 shows the sequence of processes involved in our adversarial-Q-learning scheme. Here it stands for the iteration number and s is the number of sessions. We use a greedy action in which the action corresponding to the maximum value in the Q matrix, for the state of the microswimmer, is performed; -greedy step ensures with probability g that the non-optimal action is chosen. Furthermore, we find that episodic updating of the values on the Q matrix lead to a deterioration of performance; therefore, we use continuous updating of Q. From our DNS of the 3D Navier-Stokes equation we obtain a statistically steady, homogeneous-isotropic turbulent flow in a 128 × 128 × 128 periodic domain. We introduce passive microswimmers into this flow. To define the states, we fix a coordinate triad, defined by T, ( T × ω), T⊥ as shown in fig. 9; here, T is the unit vector pointing from the microswimmer to the target, ω is the vorticity pseudo-vector, and T⊥ is defined by the conditions T⊥ • T = 0 and T⊥ • ( T × ω) = 0. This coordinate system is ill-defined if T is parallel to ω. To implement our Q-learning in 3D, we define 13 states: S = (S |ω| , S θ , S φ ) (see fig. 10); and 6 actions, A = T, -T, ( T × ω), -( T × ω), T⊥ , -T⊥ . Consequently, the Q matrix is an array of size 13 × 6.

FIG. 9. We define a Cartesian coordinate system by using the ortho-normal triad T, ( T × ω), T⊥ ; thus, all the vectorial quantities are represented in terms of this observerindependent coordinate system.