{'fraction_non_alphanumeric': 0.04173525823987115, 'fraction_numerical': 0.012814287575352037, 'mean_word_length': 4.42765314649354, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Kinetic theory of point vortices: diffusion coefficient and systematic drift

We develop a kinetic theory for point vortices in two-dimensional hydrodynamics. Using standard projection operator technics, we derive a Fokker-Planck equation describing the relaxation of a ``test'' vortex in a bath of ``field'' vortices at statistical equilibrium. The relaxation is due to the combined effect of a diffusion and a drift. The drift is shown to be responsible for the organization of point vortices at negative temperatures. A description that goes beyond the thermal bath approximation is attempted. A new kinetic equation is obtained which respects all conservation laws of the point vortex system and satisfies a H-theorem. Close to equilibrium this equation reduces to the ordinary Fokker-Planck equation.


I. INTRODUCTION
It is often useful in two-dimensional turbulence to approximate a continuous vorticity field by a cloud of point vortices. The main interest is that such a system is Hamiltonian [1] and can be studied by rather ordinary statistical mechanics. This was first considered by Onsager [2] who showed qualitatively the existence of negative temperature states at which vortices cluster. He could therefore explain the occurence of large scales vortices (or "supervortices") often observed in nature. This was a remarkable anticipation since observations were very scarce at that time. His work was pursued by Joyce & Montgomery [3] and Lundgren & Pointin [4] who introduced a mean-field approximation and obtained explicit results for the equilibrium state. They derived in particular a Maxwell-Boltzmann statistics for the distribution of point vortices at equilibrium.
Less is known concerning the relaxation towards equilibium. In fact, the evolution of the N-particle distribution function is governed by a Liouville equation but this equation contains too much information to be of practical use. One is more interested by the evolution of the one-particle distribution function P (r, t) which gives the probability that a point vortex be found in r at time t. In Ref. [5], we have described the relaxation of P (r, t) towards the Boltzmann distribution in terms of a phenomenological Fokker-Planck equation.
In this approach, the vortices have a diffusive motion due to random fluctuations and they experience in addition a systematic drift  directed along the background density gradient. Physically, the drift is the result of a polarization process and its mathematical expression can be determined with a linear response theory [5]. It is found that the drift is "attractive" at negative temperatures so the vortices cluster into macrovortices in agreement with Onsager's thermodynamical approach. At equilibrium, the drift balances the scattering and maintains a non trivial vortex distribution (the Boltzmann distribution) providing a dynamical explanation for the persistence of clustering.
In this paper, we justify our phenomenological model by deriving the Fokker-Planck equation directly from the Liouville equation, using projection operator technics [6]. These methods are standard in statistical mechanics but they are applied here for the first time to a system of vortices. We first consider the relaxation of a "test" vortex in a thermal bath in which the "field" vortices are in statistical equilibrium. In this approximation, the Fokker-Planck equation appears in its usual form with a diffusion term and a drift term. The drift coefficient is connected to the diffusion coefficient and to the temperature of the bath 1/β by an Einstein formula. The diffusion coefficient is expressed as a Kubo formula, i.e. as the integral of the velocity autocorrelation function. Using an approximation in which the vortices are advected by the equilibrium flow, we find that the autocorrelation function decays like t −2 for large times [5]. This is a slow decay but it ensures the convergence of the diffusion coefficient. We also derive non Markovian equations that keep track of memory effects. Then, we relax the thermal bath approximation and derive a generalized kinetic equation for our vortex system. This integrodifferential equation satisfies all conservation laws of the point vortex system and increases the Boltzmann entropy (H-theorem). The relaxation towards equilibrium is due to a condition of resonance. If this condition is not satisfied, the system can remain frozen in a sort of "metastable" equilibrium. By contrast, if the system is sufficiently resonant, it will converge towards the maximum entropy state described by the Boltzmann distribution. Close to equilibrium, our generalized kinetic equation reduces to the ordinary Fokker-Planck equation.
The methods developed in this paper are inspired by those introduced in plasma physics and stellar dynamics [7]. In particular, the systematic drift [5] of a point vortex is the counterpart of the dynamical friction [8] experienced by a star in a stellar system. Further analogies between 2D vortices and stellar systems are discussed in the paper and in [9][10][11][12][13][14]. Other kinetic theories of point vortices have been developed in Refs. [15,16] in a different context. A good review on point vortex dynamics is given by Newton [17].

A. The point vortex model
In a two-dimensional incompressible fluid, the velocity field u is divergenceless and can be written in terms of a stream function ψ in the form where z is a unit vector normal to the flow. The stream function is related to the vorticity ωz = ∇ ∧ u by the Poisson equation ω = −∆ψ (2) obtained by taking the curl of equation (1). The impermeability condition implies that ψ is constant on the boundary and we shall take ψ = 0 by convention. We shall consider the situation in which the velocity is created by a collection of N point vortices of equal circulation γ. In that case, the vorticity field can be expressed as a sum of δ-functions in the form where r i (t) denotes the position of point vortex i at time t. Its velocity is given by where ψ is a solution of the Poisson equation (2) with the vorticity field (3). In an unbounded fluid, one has Therefore, the velocity of a point vortex is equal to the sum of the velocities V(j → i) produced by the N − 1 other vortices, i.e. with The above dynamics can be cast in a Hamiltonian form [1]: where the coordinates (x, y) of the point vortices are canonically conjugate. These equations of motion still apply when the fluid is restrained by boundaries, in which case the Hamiltonian (9) is modified so as to allow for vortex images, and may be constructed in terms of Green's functions depending on the geometry of the domain. Since H is not explicitly time dependant, it is a constant of the motion and it represents the "potential" energy of the point vortices (we shall see later on that it also represents the kinetic energy of the flow). Therefore, point vortices behave like particles in interaction like electric charges or stars. Note, however, that the Hamiltonian (9) does not involve a "kinetic" energy of the vortices in the usual sense. This is related to the particular circumstance that a point vortex produces a velocity not an acceleration. As a result, an isolated vortex remains at rest contrary to a material particle which has a rectilinear motion due to its inertia. Point vortices form therefore a very peculiar Hamiltonian system.

B. The microcanonical approach of Onsager
The statistical mechanics of point vortices was first considered by Onsager [2] who showed the existence of negative temperatures at which point vortices cluster into "supervortices". Let us briefly recall his argumentation.
Consider a liquid enclosed by a boundary, so that the vortices are confined to an area A. Since the coordinates (x, y) of the point vortices are canonically conjugate, the phase space coincides with the configuration space and is finite: This striking property contrasts with most classical Hamiltonian systems considered in statistical mechanics which have unbounded phase spaces due to the presence of a kinetic term in the Hamiltonian. As is usual in the microcanonical description of a system of N particles, we introduce the density of state g(E) = dx 1 dy 1 ...dx N dy N δ E − H(x 1 , y 1 , ..., x N , y N ) (11) which gives the phase space volume per unit interaction energy E. The phase space volume which corresponds to energies H less than a given value E can be written It increases monotonically from zero to A N when E goes from E min to +∞. Therefore, g(E) = dΦ(E)/dE will have a maximum value at some E = E m , say, before decreasing to zero when E → +∞.
In the microcanonical ensemble, the entropy and the temperature are defined by For E > E m , S(E) is a decreasing function of energy and consequently the temperature is negative. Now, high energy states E ≫ E m are clearly those in which the vortices are crowded as close together as possible. For energies only slighlty greater than E m , the concentration will not be so dramatic but there will be a tendency for the vortices to group themselves together on a macroscopic scale and form "clusters" or "supervortices". By contrast, for E < E m , the temperature is positive and the vortices have the tendency to accumulate on the boundary of the domain in order to decrease their energy. For a system with positive and negative vortices, the negative temperature states, achieved for relatively high energies, consist of two large counter-rotating vortices physically well separated in the box. On the contrary when E → −∞, the temperature is positive and vortices of opposite circulation tend to pair off.

C. The mean field approximation
It is easy to show that the exact distribution of point vortices expressed in terms of δ-functions is solution of the Euler equation where u ex is the exact velocity field determined by equations (1) (2) and (14). This is proved as follows. Taking the derivative of (14) with respect to time, we obtain Using V i = u ex (r i (t), t), we can rewrite the foregoing equation in the form Since the velocity is divergenceless, we obtain Therefore, the Euler equation (15) with (1) (2) and (3) contains exactly the same information as the Hamiltonian system (8) (9). This description in terms of δ-functions, while being technically correct, is useless for practical purposes, because it requires the knowledge of the exact trajectories of the point vortices for an arbitrary initial condition (or the solution of the Euler equation (15)). When N is large, this task is impossibly difficult. Therefore, instead of the exact vorticity field expressed in terms of delta functions, one is more interested by functions which are smooth. For that reason, we introduce a smooth vorticity field ω (r, t) which is proportional to the average number of vortices contained in the cell (r, r + dr) at time t. This mean field description, which ignores the granularities of the sytem, requires that it is possible to divide the domain in a large number of cells in such a way that each cell is (a) large enough to contain a macroscopic number of point vortices but (b) small enough for all the particles in the cell can be assumed to possess the same average characteristic of the cell.
In this mean field approximation, the Hamiltonian (9) is changed into In writing this expression we have not taken into account the constraint j = i appearing in (9). Really, in (19) the integration extends over the point r = r ′ so that (19) contains selfenergy terms which become infinitely large for point vortices. As will soon become apparent, this meanfield approximation implies that the energy E is positive, a property which is not necessarily shared by the Hamiltonian (9). Using equation (5), adequatly generalized to account for a continuous distribution of vortices, our expression (19) for E can be rewritten Introducing explicitly the Poisson equation in (20) and integrating by parts, one has successively where u is the smooth velocity field (the second equality is obtained by a part integration with the condition ψ = 0 on the boundary). Therefore, E can be interpreted either as the potential energy of interaction between vortices (see equation (20)) or as the kinetic energy of the flow (see equation (21)).

D. The mean field equilibrium
We now wish to determine the distribution of vortices at equilibrium following a statistical mechanics approach [3]. To that purpose, using Boltzmann procedure, we divide the macrocells (r, r + dr) into a large number of microcells and enumerate the number of "microstates" which correspond to the same "macroscopic" configuration of the system. The logarithm of this number defines the entropy. In the mean field approximation, this leads to the classical formula where P (r) is the density probability that a point vortex be in the surface element centered on r. The average vorticity in r is related to this probability density by At equilibrium, the system is in the most probable macroscopic state, i.e. the state that is the most represented at the microscopic level. This optimal state is obtained by maximizing the Boltzmann entropy (22) at fixed energy (20) and vortex number N, or total circulation Writing the variational principle in the form: where β and α are Lagrange multipliers, we find that the maximum entropy state corresponds to the Boltzmann distribution [3]: with inverse temperature β. We can account for the conservation of angular momentum L = ω r 2 d 2 r (in a circular domain) and impulse P = ω yd 2 r (in a channel) by introducing appropriate Lagrange multipliers Ω and U for each of these constraints. In that case equation (26) remains valid provided that we replace the streamfunction ψ by the relative streamfunction ψ ′ = ψ + Ω 2 r 2 − Uy. This more general situation has been considered in, e.g., Ref. [18] to describe rotating or translating dipoles.
Substituting the Boltzmann relation between ω and ψ in the Poisson equation (2), we obtain a differential equation for the streamfunction: which determines the equilibrium distribution of vortices. In the case of stellar systems and electric charges, the corresponding Boltzmann-Poisson equation has the form [19]: where Φ denotes successively the gravitational and the electrostatic potential. For these systems, β ′ > 0 since the temperature is a measure of the kinetic energy. By contrast, for point vortices there is no kinetic term in the Hamiltonian (9) and the temperature can be either positive or negative. When β < 0, equation (27) is similar in structure to equation (28). The vortices tend to attract each other, like stars in a galaxy, and form "clusters" or "supervortices". The density profile determined by (27) or (28) is a decreasing function of the distance. When β > 0, equation (27) is similar in structure to equation (29). The vortices tend to repell each other, like electric charges, and accumulate at the boundary. The density profile determined by (27) or (29) is an increasing function of the distance. Therefore, the formal analogy between 2D vortices and stellar systems is intimately related to the existence of negative temperatures in 2D turbulence. However, the physical mechanism by which vortices and stellar systems achieve equilibrium is different. Whereas the organization of stars is relatively clear because of the attractive nature of gravity, the organization of point vortices at negative temperatures is much less intuitive. In the following section, we shall give a physical interpretation of this phenomenon in terms of a "systematic drift".

III. ELEMENTARY DERIVATION OF THE SYSTEMATIC DRIFT
A. Analogy with Brownian motion: the necessity of the drift We shall first show the necessity of this drift by using an analogy with Brownian theory. The starting point of this analogy is to realize that the velocity of a point vortex can be decomposed in two terms: a smoothly varying function of position and time V (r, t) and a function V V V(t) taking into account the "granularity" of the system and undergoing strong discontinuities. The total velocity of a point vortex can therefore be written: The velocity V (r, t) reflects the influence of the system as a whole and is generated by the mean vorticity ω (r, t) according to the Biot & Savart formula: The fluctuation V V V(t) arises from the difference between the exact distribution of the point vortices ω exact (r, t) and their "smoothed-out" distribution ω (r, t). It is on account of these fluctuations that the velocity of the test vortex will depart from its mean field value V . The velocity fluctuation V, of order γ d (where d ∼ n −1/2 is the inter-vortex distance), is much smaller than the average velocity V , of order nγR (where R is the domain size), but this term has a cumulative effect which gives rise to a process of diffusion. It makes sense therefore to introduce a stochastic description of the vortex motion like for colloidal suspensions in a liquid [20] or stars in globular clusters [8]. However, contrary to the ideal Brownian motion, point vortex systems have relatively long correlation times. This makes the study much more complicated than usual and the technical study of section IV is required. In order to gain some physical insight in the problem, we shall ignore this difficulty for the moment and describe the system by traditional stochastic processes.
According to equation (30), we would naively expect that the evolution of the density probability P (r, t) be governed by a diffusion equation of the form This would in fact be the case for a passive particle having no retroaction on the vortices or when the distribution of vortices is uniform like in [21] [13] [22]. However, this diffusion equation cannot be valid when the system is inhomogeneous. A first apparent reason is that equation (32) does not converge towards the Boltzmann distribution (26) when t → +∞. Another related difficulty is that equation (32) does not conserve energy. It seems therefore that a term is missing to act against the diffusion. These problems are similar to those encountered in Brownian theory or for stellar systems. They have traditionally be solved by introducing a dynamical friction in order to compensate for the effect of diffusion. The occurence of this frictional force is a manifestation of the "fluctuation-dissipation" theorem in statistical mechanics. In the present context, the dynamical friction is replaced by a systematic drift of the vortices. We must therefore rewrite the decomposition (30) in the form where ξ is the drift coefficient. In section III B, we shall give a physical justification for the existence of the drift in terms of a polarization process and in section IV we shall derive this term directly from the Liouville equation by using projection operator technics. The importance of this drift was first pointed out by Chavanis (1998) using a thermal bath approximation and a linear response theory. The drift term must be calculated by resorting to relatively elaborate technics but it is remarkable that a general relationship between ξ and D can be obtained without beeing required to analyze at any point the details of the "subdynamics". According to equation (33), the equation of motion for a point vortex can be written in the form Since the velocity V V V(t) undergoes strong discontinuities, the trajectory r(t) of the point vortex is not differentiable. Therefore, equation (34) must be viewed as a stochastic equation analogous to the Langevin equation in the ordinary Brownian theory. Let ∆t be an interval of time long compared to the fluctuation time but short at the scale on which the physical parameters change appreciably. The variation in the position ∆r of the particle during ∆t is given by where Each fluctuation V V V produces a small displacement δr but the repeated action of these fluctuations produces a net displacement of the same order as the drift ξ∇ψ. To determine the probability w[B(∆t)] that the fluctuations produce a displacement B(∆t) during the time interval ∆t, we first divide the interval (t, t + ∆t) into a succession of discrete increments in position and observe that B(∆t) is a sum of N random variables T (V i )V i where T (V) characterizes the typical duration of the velocity fluctuation V V V. This is a problem of random walks where B(∆t) represents the distance reached after N steps. For large N's, the Central Limit Theorem leads to a Gaussian transition probability: with a diffusion coefficient We now assume that the motion of a point vortex can be idealized by a Markov process, i.e. the probability at time t + ∆t depends on the probability at time t but not at earlier times. As indicated previously, this approximation is not completely correct in the case of point vortices which have long correlation times. However, relaxing this hypothesis would involve more intricate equations (see section IV) and we shall ignore this difficulty for the moment. We write therefore: where w(r − ∆r|∆r) is the probability for a point vortex located in r − ∆r to suffer an increment of position ∆r during ∆t. Expanding P (r, t + ∆t), P (r − ∆r, t) and w(r − ∆r|∆r) in the form of Taylor series, we arrive at the Fokker-Planck equation in its general form: where ∆r i ∆r j = ∆r i ∆r j w(r|∆r)d 2 (∆r).
According to equation (37), the transition probability from r to r + ∆r is given by With (43), the moments (41) (42) can be easily evaluated yielding Substituting these results in the general Fokker-Planck equation (40), we find that: We had previously introduced this equation in Ref. [5] using phenomenological arguments. The physical interpretation of each term is straightforward. The left hand side (which can be written dP/dt) is an advection term due to the smooth mean field velocity V . The right hand side can be written as the divergence of a current −∇.J and is the sum of two terms: the first term is a diffusion due to the erratic motion of the vortices caused by the fluctuations V V V; the second term accounts for the systematic drift of the vortices. At equilibrium, the drift precisely balances random scatterings and the distribution (26) is settled. More precisely, the condition that the Maxwell-Boltzmann statistics (26) satisfies equation (45) identically requires that D and ξ be related according to the relation which is a generalization of the Einstein formula to the case of point vortices. A more rigorous justification of this relation will be given in section IV where the diffusion coefficient and the drift term are calculated explicitly.

B. Systematic drift: the result of a polarization process
According to the previous discussion, the relaxation of a point vortex towards statistical equilibrium can be described by a Fokker-Planck equation involving a diffusion term −D∇P and a drift term The drift is normal to the mean field velocity V = −z∧∇ψ of the vortices and its direction, depending on the sign of β, has important physical implications. To fix the ideas, let us assume that all point vortices have positive circulation (the opposite case leads to the same conclusion). Due to the meanfield velocity, a particular point vortex rotates anticlockwise. At negative temperatures, the drift is directed to its left and the vortex is attracted to the center of the domain. On the contrary, at positive temperatures, the drift is directed to its right and the vortex is rejected against the boundary. Therefore, the effect of the drift is consistent with Onsager thermodynamical approach and it provides, in addition, a physical mechanism to understand the organization of point vortices at negative temperatures. For β = 0, the medium is homogeneous and there is no drift. Equation (47) reduces to a pure diffusion equation like in Ref. [21,13,22]. Therefore, the drift occurs only in the presence of a background shear.
In stellar systems, the relaxation of the distribution function f (r, v, t) is usually described by the Kramers-Chandrasekhar equation: which is a particular Fokker-Planck equation with a structure analogous to equation (47). In this analogy, the dynamical friction experienced by a star as a result of close encounters (Chandrasekhar [8]): is the counterpart of the systematic drift (48) experienced by a point vortex in two dimensional turbulence. The dynamical friction can be viewed as the drag exerted on a test star by the wake it induces in the field stars, like in a polarization process. We can use a similar approach to understand the origin of the drift. Let us consider a collection of N point vortices at statistical equilibrium with inverse temperature β. When β < 0, the density of these "field" vortices decreases from the center to the periphery of the domain. A "test" vortex moving through this medium locally modifies the vorticity distribution and produces a polarization cloud which amounts typically to a rotation of the surrounding vortices. This creates an excess of density behind it an a deficit of density in front of it. Therefore, the retroaction of the field vortices leads to a drift of the test vortex directed inward. We reach the opposite conclusion if the temperature is positive. When the system is homogeneous (β = 0 in a domain with no special symmetry), the polarization cloud induced by the test vortex has no effect and the drift cancels out. Therefore, a homogeneous system of point vortices remains homogeneous [21,13,22].

C. The Maximum Entropy Production Principle
In section IV, we shall derive the Fokker-Planck equation (47) directly from the Liouville equation. However, we want to show first that the general structure of this equation can be understood from relatively simple thermodynamical arguments.

Let us rewrite equation (47) in the form
where J is an unknown diffusion current. This equation conserves the circulation (24) provided that J.n = 0 on the domain boundary (n is a unit vector normal to the boundary).
The problem at hands consists in determining the expression for J. Its exact expression depends on microscopic processes and is therefore difficult to capture. However, it is easy to write down some macroscopic constraints that it must satisfy. These constraints are provided by the first and second principles of thermodynamics, namely the conservation of energy and the increase of entropy. We shall find that these constraints are very stringent and determine completely the structure of the diffusion current.
Taking the time derivative of equations (20) (22) and substituting for (51) we obtain the constraintsĖ We shall now introduce an optimization procedure known as the Maximum Entropy Production Principle (M.E.P.P.). This principle was introduced initially in the context of 2D turbulence by Robert & Sommeria [23] but its domain of applicability is very general and concerns, for example, the case of stellar systems [9,24]. This principle states that "out of equilibrium, the system evolves so as to maximize its rate of entropy productionṠ while accounting for all the constraints imposed by the dynamics, in particular the conservation of energyĖ = 0". There is no precise justification for this principle and it is important therefore to confront the M.E.P.P. with more rigorous methods, like those of section IV, to determine its domain of validity. In any case, the M.E.P.P. can be considered as a convenient tool to build relaxation equations which are mathematically well-behaved and which can serve as numerical algorithms to calculate maximum entropy states. We seek therefore the optimal diffusion current J which maximizes the rate of entropy productionṠ at fixed energy. In order to avoid the unphysical solution |J| → +∞ with J.∇ψ = 0, we impose the additional constraint where C is an upper bound which must exist but is not known. The solution of the optimization problem is where β(t) and D(r, t) are Lagrange multipliers associated with the constraints (52) and (54). When substituted in equation (51), we obtain: This equation has the same form as the Fokker-Planck equation (47). Here, the diffusion term arises from the variations of entropy δṠ and the drift term is necessary to conserve energy. Note that the Einstein formula (46) is automatically satisfied by this variational approach.
The time evolution of the inverse temperature β(t) is determined by the conservation of energy. Substituting the diffusion current (55) in the constraint (52) we find We can also check that the entropy monotonically increases during the relaxation provided that D ≥ 0. Indeed, using (53)(55) and (52) we can easily establish thaṫ At equilibrium, J = 0 and we recover the Boltzmann distribution (26).
Note that the optimal current (55) can be written J = χ∇α where α = ln ω + βγψ is a "generalized potential" which is uniform at equilibrium. Therefore, the M.E.P.P. can be viewed as a variational formulation of the linear thermodynamics of Onsager which relates the diffusion currents to the gradients of generalized potentials. However, the M.E.P.P. gives a more elegant approach to the problem and, conceding a real importance to the constraints, it is easier to implement in more complicated situations [23,9,25]. In addition, it shows that the structure of the relaxation is determined by purely thermodynamical arguments. All explicit reference to the subdynamics is encapsulated in the diffusion coefficient which is left unspecified (it appears as a Lagrange multiplier related to an unknown bound on the diffusion current). It must be therefore calculated with a more microscopic model like the one of section IV.

A. The Liouville equation
Let us consider a collection of N + 1 point vortices with identical circulation γ. We select one of these vortices, for example point vortex 0, and call it the "test vortex". The other vortices 1, ..., N will be refered to as the "field vortices". Let µ(r, r 1 , ..., r N , t) denote the N + 1 particle distribution of the system, i.e µ(r, r 1 , ..., r N , t)d 2 rd 2 r 1 ...d 2 r N represents the probability that point vortex 0 be in the cell (r, r + dr), point vortex 1 in the cell (r 1 , r 1 + dr 1 )... and point vortex N in the cell (r N , r N + dr N ) at time t. The (N + 1)-particle distribution function µ(t) satisfies the Liouville equation where V i is the velocity of vortex i produced by the other vortices according to equations (6) (7). We also introduce the one-and N-particle distribution functions defined by The physical picture that we have in mind is that the test vortex evolves in a "bath" of field vortices. Therefore, we rewrite the distribution function µ in the suggestive form where µ I reflects the effect of correlations between the test vortex and the field vortices. Physically, this term accounts for the polarization process described qualitatively in section III B. The Liouville equation (59) provides the correct starting point for the analysis of the dynamics of our vortex system. However, when N is large, this equation contains much more information than one can interpret. Consequently, what one would like to do is to describe the system in some average sense by a one-particle distribution function. In the previous sections we have indeed derived heuristically some differential equation satisfied by this distribution function on the basis of stochastic arguments. We shall now discuss the connection of such heuristic theories with a more microscopic description of the system.

B. The projection operator formalism
Our first objective is to derive some exact kinetic equations satisfied by P (r, t) and µ sys (r 1 , ..., r N , t). This can be achieved by using the projection operator formalism developed by Willis & Picard [6]. This formalism was also used by Kandrup [7] in the context of stellar dynamics to derive a generalized Landau equation describing the time evolution of the distribution function of stars in an inhomogeneous medium. We shall just recall the main steps of the theory. More details can be found in the original paper of Willis & Picard [6] and in Kandrup [7]. To have similar notations, we set x ≡ {r} and y ≡ {r 1 , ..., r N }. Then, equation (62) can be put in the form where we have written f (x, t) ≡ P (r, t) and g(y, t) ≡ µ sys (r 1 , ..., r N , t). The Liouville equation is also cast in the form where L 0 and L sys act respectively only on the variables x and y, whereas the interaction Liouvillian L ′ acts upon both x and y (the complex number i is here purely formal and has been introduced only to have the same notations as Ref. [6,7]).
Following Willis & Picard, we introduce the time-dependant projection operator: We can easily check that: We also verify that P is a projection in the sense that P 2 (t) = P (t). Applying P and 1 − P on the Liouville equation (65), we obtain the coupled equations and These equations should be compared with equations (8) (9) that appear in the quasilinear theory of 2D turbulence [26]. In the present context, equations (69) (70) describe the separation between a "macrodynamics" and a "subdynamics". In the quasilinear theory, equations (8) (9) describe the evolution of the "coarse-grained" and "fine-grained" components of the vorticity.
Introducing the Greenian we can immediately write down a formal solution of equation (70), namely: where we have assumed that initially the particles are uncorrelated so that µ I (x, y, 0) = 0. Substituting for µ I (x, y, t) from equation (72) in equation (69), we obtain The integration over y will yield an equation describing the evolution of f . Using some mathematical properties of the projection operator (66), the final result can be put in the nice symmetrical form given by Willis & Picard [6]: where the notations stand for Similarly, after integrating over x we find the equation satisfied by g [6]: C. Application to the point vortex system The previous theory is completely general and we now consider its application to a system of point vortices. Let us first rewrite the Liouville equation (59) in a form that separates the contribution of the test vortex from the contribution of the field vortices: The different operators that arise in the decomposition (65) are The mean-field velocity created by the field vortex i on the test vortex is denoted by Similarly, denotes the mean field velocity created by the test vortex on the field vortex i. Finally, the total mean field velocity experienced by the test vortex is given by where the second equality follows from the identity of the point vortices.
We are now ready to evaluate the quantities (75) (76) (77). After straightforward integration by parts, we find successively Introducing the velocity fluctuations: we can rewrite our expression for ∆ t L ′ in the form Substituting these results in equation (74), we obtain the following kinetic equation for the one-particule distribution function of a vortex system: or, alternatively, where the Greek indices refer to the components of V V V in a fixed system of coordinates. We can note that equation (93) already shares some analogies with the Fokker-Planck equation of section III. Indeed, the first term on the r.h.s. corresponds to a diffusion and the second term to a drift. For a passive particle V ν (0 → j) = 0 and the drift cancels out, as expected.

D. The thermal bath approximation
Equation (93) is an exact differential equation for P (r, t). However, this equation is not directly soluble since the unkown function µ sys ({r k }, t) is given by an equation of the form (78) depending in turn on P (r, t). We therefore have to solve the coupled system (74)-(78). This system bears exactly the same information as the initial Liouville equation (65) and, without further simplification, is untractable.
To reduce the complexity of the problem, we shall implement a "thermal bath approximation". We assume that the field vortices are in statistical equilibrium with inverse temperature β eq . Therefore, the N-particle distribution function µ sys ({r k }) can be approximated by a product of N one-particle distribution functions P eq given by the Maxwell-Boltzmann statistics (26). In other words, we make the approximation where ψ eq is solution of the Poisson equation (2) with the equilibrium vorticity ω eq = NγP eq . Substituting explicitly for the Boltzmann distribution from equation (95) in equation (93), we obtain where . eq denotes the average with respect to the equilibriun distribution P eq . Explicating the action of the Greenian, we can rewrite our equation (96) in the form where the retarded velocity V V V(j → i, t − τ ) must be viewed as an explicit function of time.
where r i (t − τ ) denotes the position at time t − τ of the i-th point vortex located in r i ≡ r i (t) at time t. The trajectories of the point vortices between t − τ and t are determined by the complicated Greenian G(t, t − τ ) defined in (71). We need therefore to solve the exact Kirchhoff-Hamilton equations of motion (8). In fact, to a good approximation, we can consider that the point vortices are purely advected by the equilibrium mean field velocity V eq . Indeed, when N → ∞, we have already indicated that the velocity fluctuation V is much smaller than the mean field velocity V eq . Therefore, we can replace the exact Greenian G by a smoother Greenian G eq which would be obtained if the point vortices were moving in the velocity field created by the equilibrium distribution function µ eq . Formally, this Greenian is constructed with the averaged Liouville operator L eq ≡ N i=0 V i eq ∂ ∂r i . In this approximation, the correlations involving two different vortex pairs vanish and the equation can be simplified considerably. Using the results of Appendix B (see in particular equation (B11)), we find where we have used P eq (r i (t − τ )) = P eq (r i ) since P eq = f (ψ eq ) is constant along a streamline and the particles are assumed to follow the streamlines in a first approximation. Since the vortices are identical we also have Noting that the integral is dominated by the divergence of the product V µ V ν when r 1 → r, we can make a "local approximation" and replace ∂ψeq ∂r ν 1 (r 1 ) and P eq (r 1 ) by their values taken in r. For the same reason, we can neglect vortex images and replace the Kernel V(0 → 1) by its singular part (7) satisfying V(0 → 1) = −V(0 → 1). With these approximations, the kinetic equation (99) takes the form: E. The Fokker-Planck equation

Unidirectionnal flow
We will now see how the previous equation can be simplified for particular equilibrium flows. We shall first consider the case of an unidirectional flow V eq = V eq (y)x produced by a vorticity distribution ω eq (y). If we restrict ourselves to solutions of the form P = P (y, t), the kinetic equation (100) becomes where we have used y(t − τ ) = y(t) = y since the point vortices follow the streamlines of the equilibrium flow. This equation can be rewritten in the form where is the velocity autocorrelation function. In Appendix C 1, it is found that: where is the local shear of the flow, equal here to the vorticity. Note that the velocity autocorrelation function decays like τ −2 for τ → +∞ [5]. Equation (102) is a non Markovian equation since the probability P (y, t) at time t depends on the probability P (y, t − τ ) at earlier times through an integration over τ . Accordingly, the present study which explicitly takes into account memory effects is more general than the stochastic model presented in section III A. However, if we implement a Markov approximation and replace P (y, t − τ ) by P (y, t), we recover the Fokker-Planck equation of section III: ∂P ∂t = ∂ ∂y D ∂P ∂y + β eq γP ∂ψ eq ∂y (106) with a diffusion coefficient and a drift term The drift coefficient is given by an Einstein relation as expected from the general considerations of section III A. The diffusion coefficient is expressed as a Kubo formula, i.e. as the integral of the velocity correlation function (see Appendix B). In Appendix C 1 it is found that (see also Ref. [5]): The reason for the logarithmic divergence is explained in Appendix A and in Ref. [13].
We now consider the case of an axisymmetrical equilibrium flow such that V eq = V eq (r)ê θ . This flow is generated by an equilibrium vorticity field ω eq (r). If we restrict ourselves to solutions of the form P = P (r, t), the kinetic equation (100) simplifies in where we have used r(t − τ ) = r(t) = r since the point vortices follow the streamlines of the equilibrium flow. Furthermore, V r(t) and V r(t−τ ) denote the radial components of V(1 → 0) at times t and t − τ in a polar system of coordinates (ê r (t),ê θ (t)) moving with the test vortex. Equation (110) can be rewritten where C(τ ) is the velocity autocorrelation function In Appendix C 2 it is found that where Σ(r) = r d dr is the local shear of the flow. If we ignore memory effects, we obtain the Fokker-Planck equation: with a diffusion coefficient and a drift term Explicitly, the diffusion coefficient has the form (see Appendix C 2 and Ref. [5]):

The general case
We can show in the general case that the relaxation of the test vortex is described by the Fokker-Planck equation with a diffusion coefficient where |Σ(r)| = 2 −Det(Σ) is the local shear of the flow and and Det(Σ) the determinant of the stress tensor Σ µν (see section V C). In the regions where the shear cancels out, our approximations clearly break up. In particular, we cannot calculate the Kubo integral by assuming that the vortices follow the streamlines of the equilibrium flow. This is because, for a local solid rotation, the vortices always remain at the same relative distance and the correlation time is infinite. In that case it is necessary to take into account the dispersion of the vortices. An alternative derivation of the diffusion coefficient can be obtained by analyzing the statistics of velocity fluctuations created by a random distribution of point vortices [13]. When the differential rotation of the vortices is neglected (which corresponds to the opposite limit of that leading to equation (120)), we obtain (see Ref. [13] and Appendix A): Clearly, a more complete study should take into account simultaneously the effect of the shear and the dispersion of the vortices to match the two formulae (120) and (121). For t → +∞, the distribution function P (r, t) of the test vortex converges towards the Maxwell-Boltzmann statistics (26). The time of relaxation corresponds typically to the time needed by the test vortex to diffuse over a distance R, the system size. Therefore t relax ∼ R 2 /D where D ∼ γ ln N is the order of magnitude of the diffusion coefficient given by equation (120). Using Γ = Nγ and introducing the dynamical time t D ∼ ω −1 ∼ R 2 /Γ, we obtain the estimate: Since the statistical description is expected to yield relevant results for large N, we conclude that the relaxation of point vortices towards the Boltzmann distribution is a very slow process. It is plausible that a more violent relaxation be at work in the system. This problem is discussed more specifically in the conclusion.
In the previous calculations, we have assumed that the distribution of the field vortices is given by the Maxwell-Boltzmann statistics (26) which corresponds to statistical equilibrium. In the case of an arbitrary background distribution P eq , the expression of the drift is where D is still given by equation (120). Since D > 0, the drift is always directed along the density gradient. The estimate of the time of relaxation is not changed in this more general situation.

V. A GENERALIZED KINETIC EQUATION
In the previous section, we have described the relaxation of a test vortex in a "bath" of field vortices at statistical equilibrium. We would like now to relax this thermal bath approximation in order to obtain a generalized kinetic equation describing the evolution of the whole system.

A. The factorization hypothesis
If the vortices are initially decorrelated then, for sufficiently short times, they will remain decorrelated. This means that the (N + 1)-particle distribution function can be factorized in a product of (N + 1) one-particle distribution functions: If we integrate the Liouville equation (59) Therefore, for sufficiently short times, the average vorticity ω satisfies the 2D Euler equation. However, at later times, the distribution function µ differs from the pure product (124) and the Euler equation does not provide a good approximation anymore. In section IV A we have determined an exact equation (93) satisfied by the one-particle distribution function at any time. This equation is not closed, however, since it involves the N-vortex distribution function µ sys . We shall close the system by assuming that µ sys can still be approximated by a product of N one-particle distribution functions: but, contrary to section IV D, the one-particle distribution function P (r k , t) is not ascribed to the equilibrium value P eq . Physically, the decomposition (62) (126) assumes that the correlations that develop between the vortices (term µ I ) are due uniquely to the polarization cloud imposed by each individual vortex. Without this polarization, the vortices would be uncorrelated (term µ sys ). In particular, this decomposition does not take into account threebody encounters which can play a crucial role in the dynamics of vortices (in particular when the system is neutral and homogeneous, see Ref. [27,22]). These high order correlations develop on longer time scales and may be neglected in a first approach. The approximation (126) introduced in equation (93) leads to a generalized kinetic equation Repeating the steps leading from equation (96) to equation (99), we obtain where P = P (r, t) and P 1 = P (r 1 , t). We also recall that between t and t − τ the trajectories of the particles are determined from the smooth velocity field created by the vorticity distribution ω = NγP (r, t). This non Markovian integrodifferential equation is similar to the equation obtained by Chavanis (2000) [26] in two-dimensional turbulence using a quasilinear theory of the Euler-Poisson system. It is shown in Appendix D that this equation rigorously conserves angular momentum in a circular domain and linear impulse in a channel (or in an infinite domain). However, under this form, it is not possible to prove the conservation of energy and the H-theorem. In our previous investigation of the problem [26], we have considered a somewhat crude approximation which amounts to neglecting memory effects in equation (128). We beleived that this approximation would not introduce any significant error but, as we shall see, we were wrong: such approximation breaks the conservation of energy.
If we assume that the decorrelation time τ is short (which does not need to be the case) and implement a strong Markov approximation, we obtain In the case of an infinite domain V(0 → 1) = −V(0 → 1) and we have the further simplification where and ξ ξ ξ = r 1 − r (we also recall that ξ ξ ξ ⊥ is the vector ξ ξ ξ rotated by + π 2 ). To arrive at equation (131) we have explicitly used the form of the Kernel (7) and to get the second equality we have used the fact that we are in two dimensions. Note that the symmetrical form of equation (130) is reminiscent of the Landau equation introduced in plasma physics and in stellar dynamics [7]. In this analogy, the position r of the vortices plays the role of the velocity v of the electric charges or stars. Therefore, we can directly infer the conservation of linear impulse P ⊥ = ωrd 2 r and angular momentum L = ωr 2 d 2 r which play respectively the role of impulse P = f vd 3 v and kinetic energy K = f v 2 2 d 3 v in plasma physics. In addition, we can prove a H-theorem for the Boltzmann entropy (22)  which is the maximum entropy state at fixed circulation, angular momentum and impulse. It is in general different from the Boltzmann distribution (26) with the relative streamfunction ψ ′ = ψ + Ω 2 r 2 − U ⊥ .r. This clearly indicates that equation (129) does not conserve the meanfield energy (20). It may happen, however, that the energy is approximately conserved if we start from an initial condition with a value of energy E 0 corresponding to β → 0, Ω → +∞ and α = βΩ/2 finite at equilibrium. In that case, equation (132) is the maximum entropy state at fixed Γ, L and E 0 (the conservation of impulse can be satisfied trivially by taking the center of vorticity as origin of our system of coordinates). However, if E differs from E 0 by a large amount, the kinetic equation (129) will not correctly describe the evolution of the system for late times. Now, if we account properly for memory effects in equation (128), we can obtain a more general equation which guaranties in addition the conservation of energy and is therefore more satisfactory. In the case of an axisymmetrical flow, it is possible to calculate the memory function appearing in equation (128) explicitly if we assume that the correlation time is smaller (but not necessarily much smaller) than the typical time on which the average vorticity changes appreciably. In this approximation, the point vortices follow, between t and t−τ , circular trajectories with angular velocity Ω(r, t) = V θ (r, t)/r and equation (128) simplifies in (see Appendix E): where Ω = Ω(r, t), Ω 1 = Ω(r 1 , t) and r > (resp. r < ) is the biggest (resp. smallest) of r and r 1 . The angular velocity is related to the vorticity by Similarly, in the case of a unidirectional flow, we can assume that, between t and t − τ , the point vortices follow linear trajectories with velocity V = V (y, t)e x . This leads to the kinetic equation (see Appendix E): where V = V (y, t) and V 1 = V (y 1 , t). The function E 1 (x) is the exponential integral and L an upper cut-off necessary in that case (see Appendix E). We also recall that the average velocity is related to the vorticity by We can remarkably propose an approximation of the general kinetic equation (128) which encompasses both the axisymmetric form (133) and the unidirectional form (135). Memory effects are not neglected, unlike in Eq. (130), but they are simplified in a way which preserves all the conservation laws of the system (as discussed below). We propose the generalized kinetic equation: and For specific applications, it may be necessary to introduce a shielding of the form (E39) in the interaction between vortices. This shielding arises naturally in geophysics in the "quasigeostrophic approximation". In that case the tensor K µν is replaced by (1/L)K 1 (ξ/L)ξ µ ⊥ ξ ν ⊥ /ξ where K 1 is the modified Bessel function of first order and L is an upper-cut off (called the Rossby radius in geophysics) which plays the same role as the Debye length in plasma physics.
Equation (137) is not exact, in a strict sense, yet it satisfies all the conservations laws of the vortex system as well as a H-theorem (section V B). This is very gratifying and this can have important practical applications. It is remarkable that we can write down an approximate kinetic equation in the general case without being required to specify the trajectory of the point vortices between t and t − τ . In fact, to arrive at equation (137) we have made implicitly two approximations: (i) we have assumed that the vorticity field does not change dramatically when we follow the vortices in their motion between t and t − τ . (ii) After the time integration has been effected, the non universal function V ν (1 → 0, t − τ ) gives rise to a logarithmic term which has been replaced by 1 in the subsequent calculations. This term produces a sub-logarithmic correction that is flow-dependant and that has been neglected. It is on account of this weak dependance that a general kinetic equation can be obtained. For axisymmetrical or unidirectionnal flows, equation (137) reduces to equations (133) and (135) with ln(|r 1 − r|/(r 1 + r)) and K 0 (|y − y 1 |/L) instead of ln(1 − (r < /r > ) 2 ) and E 1 (2|y − y 1 |/L). This is not a too severe discrepency (these functions have a similar logarithmic behaviour) so our approximations are reasonable.
From equations (133) (135) and (137), it is clear that the relaxation towards equilibrium is due to a phenomenon of resonance. Only the points r 1 satisfying the condition ξ ξ ξ.v = 0 with r 1 = r will contribute to the diffusion current in r. In the axisymmetrical case, this condition of resonance reduces to Ω(r 1 ) = Ω(r) and in the unidirectional case to V (y 1 ) = V (y). These conditions of resonance had never been noticed previously. Further work on the subject will have to make these criteria more precise by computing explicitly "resonance lines" in two-dimensional real flows.

B. Conservation laws and H-theorem
We now derive the conservation laws and the H-theorem satisfied by equation (137). First of all, the conservation of the circulation is straightforward since equation (137) can be written in the form of a continuity equation (D1). To prove the conservation of angular momentum, we start from equation (D4), substitute for (137), permut the dummy variables r and r 1 and add the resulting expressions. This yieldṡ But, from equation (138), we immediately verify that This proves the conservation of angular momentum. We can prove the conservation of linear impulse in a similar manner. Starting from equation (D16), substituting for (137), permuting the dummy variables r and r 1 and adding the two resulting expressions yieldṡ P = 0. For the conservation of energy, we start from equation (52) and follow the by-now familiar procedure. This yieldṡ Considering the form of the tensor (138), we have When substituted in equation (141), we see that the occurence of the delta function in the kinetic equation impliesĖ = 0. Finally, for the rate of entropy production we have, according to equations (53) and (137): Permutting the dummy variables r and r 1 and adding the resulting expression to equation (143), we geṫ Now, for any vector, A µ K µν A ν = (Aξ ξ ξ ⊥ ) 2 /ξ 2 ≥ 0. This proves a H-theorem (Ṡ ≥ 0) for our kinetic equation (137). It should be emphasized that the conservation laws and the Htheorem result from the symmetry of the kinetic equation (and the condition of resonance) and not from formal Lagrange multipliers like in the thermodynamical approach of section III C. In addition the H-theorem is proved by our approach instead of beeing postulated. This is more satisfying on physical grounds.
It remains for one to show that the Boltzmann distribution is a stationary solution of equation (137). Noting that we have successively where we have used equations (140) and (142). When substituted in equation (137), we find that the r.h.s. cancels out due to the delta function. The advective term is also zero since P = f (ψ ′ ). Therefore, the distribution (145) is a stationary solution of equation (137). Note, however, that this is not the only solution. Any stationary solution satisfying ξ ξ ξ.v = 0 for any couple of points r, r 1 (with r = r 1 ) is a solution of equation (137). Physically, this implies that the system needs sufficiently strong resonances to relax towards the maximum entropy state. If this is not realized it will be frozen in a sort of "metastable" equilibrium. This may explain why the maximum entropy state is not always reached in two-dimensional turbulence. For example, for a unidirectional flow with ω positive or negative everywhere, the velocity field is monotonous (see equation (136)) and the condition of resonance cannot be satisfied. The evolution of the system will require non trivial correlations between point vortices that are not taken into account in the present approch. One would need to replace the factorization hypothesis (126) by a product of two-point or three-point correlations functions. However, it is plausible that these correlations develop on a very long time scale so it remains a matter of debate to decide whether they really are relevant for the dynamics or not.
In the context of 2D turbulence described by the Euler-Poisson system, the quasilinear theory developed by  yields instead of equation (137): where K µν is defined by (138) and ω = ω(r, t), ω ′ = ω(r ′ , t). In addition to the previous conservation laws, this equation guaranties that the coarse-grained vorticity ω remains bounded by the maximum value of the initial distribution, i.e. ω ≤ σ 0 . This equation satisfies a H-theorem for the Fermi-Dirac entropy introduced by Miller-Robert-Sommeria at equilibrium [28,29]. Our approach provides therefore another way of justifying their results from a dynamical point of view. Equation (148) is written for a single level of vorticity σ 0 , but it is possible to extend the quasilinear theory to an arbitrary distribution of vorticity levels (in preparation). Our equations should provide therefore an interesting and useful parametrization of the 2D Euler equation. It should be recalled in that respect that the usual turbulent diffusion ν∆ω introduced ad hoc in the r.h.s. of the 2D Euler equation in order to smooth out the small scales and prevent numerical instabilities breaks the conservation laws of the inviscid dynamics. This is not the case for our equation (148): not only it smoothes out the unresolved scales (as exemplified by the existence of a H-theorem) but it satisfies all the conservation laws of the inviscid dynamics and respects the invariance properties of the Euler equation (invariance by translation and rotation of the coordinates, Galilean invariance and invariance by rotation of the referential). In addition there is no free parameter in our theory except the coarse-graining mesh ǫ (or resolution scale) which depends on the situation contemplated. Different attempts had been made previously to obtain an equation satisfying all these requirements, but only partial results were obtained [23,9,25,5,26].

C. Connexion with the Fokker-Planck equation
Equation (137) can be considered as our final result but we wish to show that a direct connection with the Fokker-Planck equation of section IV can be found. Introducing a diffusion tensor and a drift term equation (137) can be rewritten in the more illuminating form: similar to the general Fokker-Planck equation (40). Note, however, that equation (151) is an integrodifferential equation since the density probability P (r, t) in r at time t depends on the value of the whole distribution of probability P (r 1 , t) at the same time by an integration over r 1 . By contrast, the Fokker-Planck equation (47) is a differential equation. The usual way to transform an integrodifferential equation into a differential equation is to make a guess for the function P (r 1 ) appearing under the integral sign and refine the guess by successive iterations. In practice we simply make one sensible guess. Therefore, if we are close to equilibrium, it seems natural to replace the function P 1 appearing in the integrals by the Boltzmann distribution This corresponds to the "thermal bath approximation" of section IV D: the vortices have not yet relaxed completely, but when we focus on the relaxation of a given point vortex (described by P ) we can consider, in a first approximation, that the rest of the system (described by P 1 ) is at equilibrium. Within this approximation, the diffusion coefficient and the drift simplify in: where we have made the local approximation. If we assume that the correlation time is short, i.e. if we replace ξ 2 δ(ξ ξ ξ.v) by τ /π 2 (compare equations (137) and (130)), we obtain In that case, equation (151) reduces to the Fokker-Planck equation found in section IV: If instead of the Boltzmann distribution (152) we use the Gaussian distribution (132), we get This equation is closely related to the Kramers-Chandrasekhar equation (49) since the drift and the friction are linear in r and v respectively. The diffusion coefficient (156) was previously obtained by [9,30,26] using phenomenological arguments. However, in these studies the correlation time τ was not specified. A first determination of τ was obtained in Ref. [5] but it was restricted to axisymmetrical or unidirectional flows. Using equation (154), we can determine the expression of the diffusion coefficient and the correlation time in the general case. Expanding the velocity difference v = V 1 − V in a Taylor series in ξ ξ ξ = r 1 − r, we obtain to first order in the expansion where is the stress tensor. It satisfies the property of symmetry Σ µν = Σ νµ . Since the flow is divergenceless, we also have Σ xx + Σ yy = 0. This suggest to introduce the notations a = Σ xx = −Σ yy and b = Σ xy = Σ yx . In terms of the stress tensor (160), the diffusion tensor (154) can be rewritten This integral can be performed easily by working in a basis where the tensor Σ µν is antidiagonal. To that purpose, we seek a tensor T such that Σ =T Σ ′ T where Σ ′ is antidiagonal and, by definition,T µν = T νµ . We also impose that T is unitary so thatT = T −1 . Then, if we denote by (M 11 , M 12 , M 21 , M 22 ) the components of a 2 × 2 matrix, we find Σ ′ = (0, b ′ , b ′ , 0) and T = (α, β, −β, α) with α 2 + β 2 = 1 and b ′ = b/(α 2 − β 2 ) = −a/2αβ. From the above results, it is also clear that b ′2 = −Det(Σ) = (a 2 + b 2 ), where Det(Σ) stands for the determinant of the matrix Σ. Now, introducing a new system of coordinates such that ξ ′ µ = T µν ξ ν , or alternatively ξ µ = T νµ ξ ′ ν , we easily check that the Jacobian of the transformation ξ ξ ξ → ξ ξ ξ ′ , i.e. the determinant of T , is equal to one. Under these circumstances, the diffusion tensor (161) can be written D =T D ′ T with where we have set |Σ(r)| = 2 −Det(Σ). Physically, this quantity represents the local shear of the flow. The components of the tensor (162) can now be determined easily. First of all, Setting ξ ′ 1 = ξ cos θ and ξ ′ 2 = ξ sin θ where ξ = ξ ′ = |r 1 − r|, we get or, equivalently, dθ sin θδ(cos θ).
As explained in Appendix A, we regularize the logarithmic divergence by introducing appropriate cut-offs at small and large scales. With the change of variables t = cos θ, we obtain By the same arguments, we find that D ′22 = D ′11 = D. On the other hands, it follows for reasons of antisymmetry by the transformation ξ ′ 1 → −ξ ′ 1 that D ′12 = D ′21 = 0. Therefore D ′µν = Dδ µν is diagonal in the basis where Σ is anti-diagonal. This remains true in any basis since D µν = T λµ D ′λσ T σν = D(T T ) µν = Dδ µν . Therefore, close to equilibrium, the diffusion is isotropic and the general expression of the diffusion coefficient is Comparing with equation (156) we find that the correlation time is given by The quantity |Σ(r)| = 2 −Det(Σ) plays a fundamental role in the theory. Clearly, this expression is invariant by a change of referential. For a unidirectional flow or an axisymmetrical flow, we recover the results of section IV E and of Ref. [5].

VI. CONCLUSION
In this paper, we have provided a systematic derivation of the kinetic equations of point vortices, applying for the first time the powerful projection operator technics to this problem. We have described how a cloud of point vortices relaxes toward the mean field statistical equilibrium, leading to a clustering into large coherent vortices. In the first part of the paper, we have focused on the relaxation of a "test" vortex in a cloud of background vortices at statistical equilibrium. The cloud of "field" vortices plays the role of a thermal bath like in other problems of statistical physics. We have shown that the test vortex undergoes a usual diffusion effect due to random fluctuations and that it also experiences a systematic drift. This drift, due to a polarization of the background vortices by the test vortex, balances the effect of diffusion at equilibrium, providing a dynamical explanation for the persistence of clustering. The drift was previously derived with a linear response theory , but the diffusion was heuristically introduced by adding a white noise effect. The present derivation systematically derives the two effects, diffusion and drift, from the same formalism, and is therefore more satisfactory. The diffusion derived here turns out to be influenced by long time correlations, so it is more complex than the usual white noise effect. When memory effects are ignored, we obtain a Fokker-Planck equation for the evolution of the one-particle distribution function. This Fokker-Planck equation can also be derived from a phenomenological Maximum Entropy Production Principle [23]. This shows that the structure of this equation is influenced more by thermodynamics (the first and second principles) than by the precise microscopic model. However, our systematic procedure starting directly from the Liouville equation provides a justification for this thermodynamical approach and specifies its range a validity. It also allows to determine explicitly the value of the diffusion coefficient which was left unspecified by the Maximum Entropy Production Principle. All these results could be tested numerically by introducing a test vortex in "sea" of vortices at statistical equilibrium and by solving the Kirchhoff-Hamilton equations of motion.
In the second part of the paper, we have attempted to describe the evolution of the whole system of vortices far from equilibrium. We have obtained a new kinetic equation (128) which incorporates a delocalization in space and time. This is therefore a non Markovian integrodifferential equation. A similar equation also occurs in the quasilinear theory of the 2D Euler equation . Within some approximation, it is possible to carry out the time integration explicitly and this yields a simpler equation (137) which only conserves a delocalization in space. This equation respects all the conservation laws of the point vortex system and satisfies a H-theorem. The relaxation is due to a condition of resonance between distant vortices. If the system is sufficiently "resonant", it will reach a maximum entropy state described by the Boltzmann distribution. However, if there are not enough resonances the evolution may stop on metastable state. Only non trivial correlations between vortices (for example three body collisions) can unfreeze the system and induce further evolution. These correlations are not taken into account in the present analysis although the projection operator formalism might still apply. It would be necessary to modify the factorization hypothesis (126) so as to account for two-body or three-body correlation functions. It is possible, however, that these correlations develop on a much longer time scale so it is not yet clear whether they are physically relevant. In any case, the approximations made in the present paper are a first step towards a rational kinetic theory of point vortices.
It is also possible that a system of point vortices undergoes a form of "violent relaxation" in its early stage. For short time scales, the correlations between point vortices have not yet developed and the average vorticity is solution of the 2D Euler equation. When the initial condition is far from equilibrium, it is well-known that the 2D Euler equation develops a complicated mixing process leading to the formation of an organized state (on a coarsegrained scale). This relaxation is quite rapid, of the order of the dynamical time t D , and the resulting equilibrium state is predicted to be a complicated superposition of Fermi-Dirac distributions respecting all the constraints of the Euler equation [28,29]. On longer time scales, of the order of N ln N t D , the correlations between point vortices develop and the system undergoes another form of relaxation, much slower. This relaxation is towards the Boltzmann distribution derived by [3] Joyce & Montgomery and [4] Lundgren & Pointin which is the true equilibrium state for a system of point vortices. The first type of relaxation, by vorticity mixing, was described in Ref. [26] using a quasilinear theory of the 2D Euler equation. The exclusion principle leading to the Fermi-Dirac statistics was explicitly shown as well as a H theorem. The second type of relaxation, due to discrete interactions between point vortices, was the object of the present paper.
It is noteworthy that a similar distinction occurs in the context of stellar systems [9]. Indeed, the relaxation of stars is a two stage process. For short time scales ∼ t D , the encounters between stars can be neglected and the distribution function is solution of the Vlasov equation (analogous to the 2D Euler equation). If the system is initially far from mechanical equilibrium, it will experience a "violent relaxation" towards a virialized state. This equilibrium is predicted to be a superposition of Fermi-Dirac statistics [31,32], like for the 2D Euler equation. Then, on a longer time scale, of the order of N ln N t D , the encounters between stars cannot be ignored anymore and will deviate the stars from their unperturbed trajectories. This collisional relaxation is usually described by a Landau or Fokker-Planck equation that converges towards the Maxwell-Boltzmann distribution at equilibrium. On even longer time scales, three-body encounters leading to the formation of binaries induce a slow evolution of the system [33]. Three body encounters (involving vortices of different sign) are also relevant in two dimensional turbulence and lead to the formation of vortex "pairs" [22]. The analogy with "binary stars" is interesting to note. It would be important to test these ideas with numerical simulations of point vortices or stars. The situation is difficult in the stellar context because a maximum entropy state does not always exist. Indeed, the system can collapse and overheat: this is the so-called "gravothermal catastrophe" [34]. This problem does not occur for point vortices and it should be possible to evidence the two (or more) successive equilibria more properly. An advantage of point vortices with respect to stars is the lower dimensionality of space (D = 2 instead of D = 3, or D = 6 in phase space) that should make numerical simulations easier. In this Appendix, we study the statistics of velocity fluctuations produced by an inhomogeneous distribution of point vortices. This study extends the calculations of [35,36,13] for a uniform medium and provides a simple framework to understand the logarithmic divergence of the diffusion coefficient. Let us consider a collection of N point vortices randomly distributed in a disk of radius R with an average density n(r). The velocity V occuring at a given location r of the flow is the sum of the velocities Φ Φ Φ i (i = 1, ..., N) produced by the N vortices: Following a procedure similar to that adopted in Ref. [13], the velocity distribution can be expressed as Here, P (r 1 ) denotes the probability of occurence of a point vortex in r 1 and by definition Introducing explicitly the vortex density n(r 1 ) = NP (r 1 ), we have we can rewrite our expression for A N (ρ ρ ρ) in the form In the limit of large N, R with fixed n(r), we can approximate the foregoing expression by Separating the real and imaginary parts of C(ρ ρ ρ), we obtain In the first integral, we find it convenient to introduce the relative separation ξ ξ ξ = r 1 − r in terms of which In Ref. [13], it was found that an important contribution to the velocity fluctuations comes from the nearest neigbour. This justifies to make the local approximation n(r + ξ ξ ξ) ≃ n(r) in equation (A12). In this approximation This quantity is closely related to the function C(ρ ρ ρ) evaluated in Ref. [13] for a distribution of vortices with uniform density n. The only difference is the presence of the local density n(r) in place of n. Therefore, we can infer directly that In fact, the local approximation is only marginally valid because, as discussed in Ref. [13], the contribution of the nearest neighbor is precisely of the same order of magnitude as the contribution of the rest of the system. This results in a logarithmic divergence in (A15) due to the weak collective behaviour of the system. In the second integral appearing in (A11), the contribution from proximate vortices vanishes by symmetry. As a result, the integral is dominated by large values of |r 1 | or, equivalently, by small values of |Φ Φ Φ|. We can therefore make the approximation sin(ρ ρ ρΦ Φ Φ) ≃ ρ ρ ρΦ Φ Φ and write In the integral, we recognize the mean-field velocity created in r by the average distribution of vortices: Hence, Substituting the explicit expression for C(ρ ρ ρ) in equation (A9), we obtain Therefore, the velocity distribution in r can be written quite generally This is the same distribution as for a uniform distribution of vortices except that the constant density n has been replaced by the local density n(r) and that the distribution is for the fluctuating velocity V V V = V − V (r). Repeating the calculations of Ref. [13], the velocity p.d.f is explicitly given by: This distribution lies at the frontier between Gaussian and Lévy laws: the core of the distribution is Gaussian while the tail decays algebraically like for Lévy laws. This is because the variance of the individual velocities (A2) diverges logarithmically so the central limit theorem is only marginally applicable. For that reason we have proposed to call this distribution the "marginal Gaussian distribution" [14]. In the strict mathematical limit N → +∞, the transition between the two regimes is rejected to infinity and the velocity p.d.f. is purely Gaussian [36]. However, the convergence towards this Gaussian distribution is so slow that in practical applications it is never reached: the algebraic tail always remain [21]. According to equations (A1) and (A2), the variance of the velocity can be written where we have made the local approximation in the second equality. This quantity diverges logarithmically at both small and large vortex separations. The divergence at small separations is a failure of our model which ignores correlations between vortices. In fact, when two vortices approach each other they can form a pair, as discussed in Ref. [21], and our mean field theory clearly breaks down. We shall account heuristically for this failure by introducing a cut-off at some minimum distance d pair ∼ (πn ln N) −1/2 [13]. The divergence at large separations is due to the unshielded nature of the interaction potential. It is therefore natural to cut the integral at R, the typical size of the system. With this regularization, we obtain Since the divergence is weak (logarithmic), the result does not depend crucially on the precise value of the cutoffs. The same expression for the variance can also be obtained from the formula if we introduce a cut-off at large velocities and use equations (A21) (A22) (A23) [to sufficient accuracy, we just need considering equation (A21)].
The diffusion coefficient can be expressed in terms of the variance (A25) and the typical correlation time τ by [13]: Using equation (A25) and the relation ω = nγ, we obtain in agreement with our result (156). The correlation time is more difficult to evaluate. If we ignore the differential rotation of the vortices, the calculations of Ref. [13] are directly applicable and lead to the expression This is the typical time needed by a vortex to cross the interparticle distance d ∼ 1/ n(r) with the velocity V 2 . With this approximation, the diffusion coefficient is given by equation (121). On the other hand, if we consider that the vortices are transported by the equilibrium flow and evaluate the diffusion coefficient with the Kubo formula (see Appendix B), we find that the correlation time is related to the local shear by equation (168), i.e.: Physically, it corresponds to the time needed by two vortices with relative velocity Σd to be stretched by the shear on a distance ∼ d. This approximation breaks up, however, when the shear is weak. In that case τ is given by equation (A29) obtained when only the dispersion of the vortices is considered. Clearly, a general formula should take into account simultaneously the effect of the shear and the dispersion of the particles.

APPENDIX B: THE KUBO FORMULA
Let us consider the diffusion of a test vortex in a "sea" of field vortices described by the equilibrium distribution The general form of the diffusion coefficient, as defined by equation (40), writes Now, the net displacement of the test vortex produced by the fluctuations of the velocity between t and t + ∆t is given by Substituting explicitly for ∆r from equation (B3) in equation (B2), we have where . eq denotes the average with respect to the equilibrium distribution (B1). Since the correlation function appearing in the integral only depends on the time difference |t ′′ − t ′ |, we also have or, alternatively, With the change of variables τ = t ′ − t ′′ , we obtain successively If the correlation function V µ (t)V ν (t − τ ) eq decays more rapidly than τ −1 , we can take the limit ∆t → +∞ to finally obtain This is the Kubo formula for our problem. Remembering that V V V denotes the fluctuation of the total velocity: we find that Now, the first quantity in bracket can be written explicitly For large N, we can make the approximation N(N − 1) ≃ N 2 . Using (B10) and (B11) the correlation function can be put in the form Since the integral is dominated by interactions involving relatively close vortices, we can make the local approximation: The expression for the diffusion coefficient then becomes For sake of brevity, we shall denote the velocity correlation function by Therefore, the Kubo formula takes the form More generally, we have These quantities are calculated explicitly in Appendix C in the case of simple flows. For ∆t → 0 (ballistic motion): and for ∆t → +∞ (diffusive motion): ln N |Σ(y)| P eq (y)∆t. (C11)

Axisymmetrical flow
In an axisymmetrical flow, the trajectory of a fluid particle takes the simple form: As indicated in section IV E 2, we are particularly interested by the r(t)r(t − τ ) component (112) of the correlation function. Let us introduce the separation δr ≡ r 1 − r between the field vortex 1 and the test vortex. In the local approximation, δr can be considered as a small quantity. Therefore we can write: With these notations, the correlation function (112) can be rewritten: Now, and In the local approximation, we can expand the last term in equation (C18) in a Taylor series in r 1 − r. This yields where Σ(r) is the local shear of the flow (114). Substituting equations (C17) and (C19) in equation (C16), we get This integral is similar to (C6), so we directly obtain equations (113) and (118).

APPENDIX D: CONSERVATION LAWS SATISFIED BY THE GENERALIZED KINETIC EQUATION
In this section, we prove some general properties satisfied by equation (128). Note first that it can be written where is the diffusion current. It is clear at first sight that equation (D1) conserves the total circulation Γ = ω d 2 r provided that J.n = 0 on the boundary. We now prove the conservation of other integral constraints depending on the domain shape.
Permuting the dummy variables r and r 1 , we geṫ Adding these two quantities, we arrive at the final expressioṅ Now, the term in brackets vanishes as shown by the following argument. Consider two point vortices in a circular (or infinite) domain. Their angular momentum is and it is conserved. This implies: We can also prove this result by a direct calculation. In an unbounded domain V(0 → 1) = −V(1 → 0) and consequently where we have used equation (7) to get the last equality. In a circular domain, the velocity V(1 → 0) is given by equation (7) plus a term V b (1 → 0) which can be determined with the method of "images". If R denotes the domain radius, we find Therefore, Noting that (ẑ ∧ r 1 ).r = −(ẑ ∧ r).r 1 and that we finally conclude that r.V(1 → 0) + r 1 .V(0 → 1) = 0.
From this identity and from equation (D7), it results that the kinetic equation (128) conserves the angular momentum in a disk and in an infinite domain, i.e.L = 0.
(ii) In an infinite domain or in a channel, the linear impulse P = r ∧ ω ẑd 2 r (D15) must be conserved (in a channel extending in the x direction, only the component P x of the linear impulse must be conserved). Taking the time derivative of equation (D15), substituting for (D1) and remembering that P is conserved by the advective term, we geṫ Substituting explicitly for the diffusion current (D2) in equation (D16), permutting the dummy variables r and r 1 , and taking the half-sum of the resulting expressions we finally obtainṖ Now, we can use the same argument as before to show that the term in brackets vanishes. Let us consider two point vortices in a channel (or in an infinite domain). Their linear impulse is and it is conserved. This implies We can also prove this result by a direct calculation. Equation (D19) is obvious in an unbounded domain since V(0 → 1) = −V(1 → 0). In a channel extending in the x direction, we need to show thatṖ x = 0, i.e. V y (1 → 0) + V y (0 → 1) = 0. Now, the velocity V y (1 → 0) is given by equation (7) plus a term V b (1 → 0) y which can be determined with the method of "images". If a denotes the width of the channel, we find that Under this form, it is clear that V b (1 → 0) y is antisymmetric under the exchange of 1 and 0, so that finally V y (1 → 0) + V y (0 → 1) = 0. (D21) From this identity and from equation (D17), it results that the kinetic equation (128) conserves the linear impulse in a channel or in an infinite domain, i.e.Ṗ = 0.

APPENDIX E: CALCULATION OF THE MEMORY FUNCTION
In this Appendix, we calculate the memory function that occurs in equation (128).
We find similarly that V r 1 (t−τ ) (1 → 0) = r r 1 V r(t−τ ) (1 → 0). Our previous assumptions also imply that P (r(t − τ ), t − τ ) ≃ P (r, t) between t and t − τ . In words, this means that the correlation time is smaller than the time scale on which the average vorticity changes appreciably. We do not assume that it is much smaller, so this approximation is not over restrictive. In that case, the diffusion current becomes where the time integral has been extended to +∞. We now need to evaluate the memory function Introducing the notations φ = θ 1 − θ and λ = 2rr 1 r 2 1 + r 2 < 1, It remains for us to evaluate the series that appears in the last expression of the memory function. Using the identities: and the definition (E8) of λ, we find that a 0 < ∞ and, for n > 0: a n = − 1 n (r 2 1 + r 2 ) − |r 2 1 − r 2 | 2rr 1 where r > (resp. r < ) is the biggest (resp. smallest) of r and r 1 . Therefore, the value of the series is This leads to the kinetic equation (133).

Unidirectional flow
If we assume that P = P (y, t), then equation (128) simplifies in Assuming that between t and t − τ the vortices follow linear trajectories with velocity V = V (y, t)e x , we have y(t − τ ) = y and x(t − τ ) = x − V (y, t)τ . Therefore, the function V y (1 → 0) at times t and t − τ takes explicitly the form and where we have introduced the notation ∆V = V (y 1 , t) − V (y, t).
We also assume that the correlation time is smaller than the time scale over which the vorticity changes appreciably. Then, P (y(t − τ ), t − τ ) ≃ P (y, t) and the diffusion current becomes J y = −N t 0 dτ dx 1 dy 1 V y (1 → 0) t V y (1 → 0) t−τ P 1 ∂P ∂y − P ∂P 1 ∂y 1 .
We now need to calculate the memory function W (ξ ξ ξ) = ln ξ → W (ξ ξ ξ) = −K 0 (ξ/L), where K 0 is the modified Bessel function of order zero and L is a length scale of the order of the system size. For small separations K 0 (z) ∼ − ln z and for large separations K 0 (z) ∼ π 2z e −z . This modification amounts to replacing the Poisson equation (2) by an equation of the form Equation (E39) is precisely what is obtained in geophysics in the "quasigeostrophic approximation". The deformation of the fluid surface introduces a shielding of the interaction between vortices on a length ∼ L, called the Rossby radius. Obviously, the Rossby radius plays the same role as the Debye length in plasma physics.
With this prescription we find that where k L = 1/L. Substituting in equation (E35), we get The integral on k y can be carried out easily, leaving the result Now, setting t 2 = 4Y 2 (k 2 x + k 2 L ), we finally obtain where is the exponential integral. For x → 0, we have E 1 (x) = −C − ln x where C = 0.57721... is Euler's constant. In conclusion, the diffusion current in the unidirectional case takes the form and it leads to the kinetic equation (135).