P Degond
email: degond@mip.ups-tlse.fr

S Motsch
email: motsch@mip.ups-tlse.fr

Large scale dynamics of the Persistent Turning Walker model of fish behavior

Keywords: Individual based model, Fish behavior, Persistent Turning Walker model, Ornstein-Uhlenbeck process, kinetic Fokker-Planck equation, asymptotic analysis, diffusion approximation AMS Subject classification: 35Q80, 35K99, 60J70, 82C31, 82C41, 82C70, 82C80, 92D50 1

This paper considers a new model of individual displacement, based on fish motion, the so-called Persistent Turning Walker (PTW) model, which involves an Ornstein-Uhlenbeck process on the curvature of the particle trajectory. The goal is to show that its large time and space scale dynamics is of diffusive type, and to provide an analytic expression of the diffusion coefficient. Two methods are investigated. In the first one, we compute the large time asymptotics of the variance of the individual stochastic trajectories. The second method is based on a diffusion approximation of the kinetic formulation of these stochastic trajectories. The kinetic model is a Fokker-Planck type equation posed in an extended phase-space involving the curvature among the kinetic variables. We show that both methods lead to the same value of the diffusion constant. We present some numerical simulations to illustrate the theoretical results.

Introduction

This paper considers a new model of individual displacement, the so-called 'Persistent Turning Walker' (PTW) model, which has recently been introduced to describe fish behavior [START_REF] Gautrais | Analyzing fish movement as a persistent turning walker[END_REF]. The fish evolves with a velocity of constant magnitude and its trajectory is subject to random turns (i.e. random changes of curvature) on the one hand and to curvature relaxation to zero on the other hand. The random changes of curvature can be interpreted as a way for the fish to explore its surroundings while relaxation to zero curvature just expresses that the fish cannot sustain too strongly curved trajectories and when the curvature becomes too large, it gets back to a straight line trajectory. The combination of these two antagonist behaviors gives rise to an Ornstein-Uhlenbeck process on the curvature. The curvature is the time derivative of the director of the velocity, while the velocity itself is the time derivative of position. The PTW process collects all these considerations into a system of stochastic differential equations.

This model is, to the knowledge of the authors, original, and has appeared for the first time in the works by Gautrais, Theraulaz, and coworkers [START_REF] Gautrais | Analyzing fish movement as a persistent turning walker[END_REF]. The present paper considers the large time and space scale dynamics of a two-dimensional particle subject to this PTW process. It rigorously shows (in the mathematical sense) that, at large scales, the dynamics of the particle can be described by a diffusion process and it provides a formula for the diffusion coefficient. To prove this result, two methods are considered.

In the first method, the stochastic differential system itself is considered and the variance of the position is shown to behave, at large times, like a linear function of time. The diffusion coefficient is identified as the slope of this linear function. Because the curvature and the velocity angle can be explicitly computed, an explicit formula for the diffusion coefficient, involving some special functions, can be obtained.

The second method considers the forward Kolmogorov equation of the stochastic process. This equation gives the evolution of the probability distribution function of the particle in the extended phase space (position, velocity angle, curvature) as a function of time. It is a Fokker-Planck type equation. The passage from the microscopic to the macroscopic scales relies on a rescaling of the Kolmogorov equation. This rescaling depends on a small parameter ε ≪ 1, which describes the ratio of the typical microscopic to macroscopic space units. After this rescaling, the problem has the typical form of the diffusion approximation of a kinetic problem (see references below). The goal is then to study the behaviour of the solution as ε → 0. It is shown that the solution converges to some 'thermodynamical equilibrium' which is a Gaussian distribution of the curvature and a uniform distribution of the velocity angle. The equilibrium depends parametrically on the density which satisfies a spatial diffusion equation.

Finally, the connection between the two methods is made by showing that the diffusion tensor in the second approach can be represented by a formula involving the solution of the stochastic differential equation of the first approach. Additionally, this representation leads to explicit computations which show that the two formulas for the diffusion coefficient actually coincide. This seemingly innocuous result is actually quite powerful. Indeed, the diffusion approximation method leads to a non-explicit expression of the dif-fusion coefficient, involving the moments of a particular solution of a stationary equation involving the leading order operator of the Fokker-Planck equation. That this non-explicit formula is equivalent to the explicit formula given by the stochastic trajectory method is by far not obvious. In this respect, the stochastic trajectory method is more powerful than the diffusion approximation approach, because it directly leads to the most simple expression of the diffusion constant.

A third route could have been taken and has been dismissed. This third method would actually use the stochastic differential equation itself to perform the diffusion approximation in the forward Kolmogorov equation. We have preferred to use partial differential equation techniques. One reason for this choice is that these techniques can be more easily extended to more complex situations. One typical example of these more complex situations are the nonlinear systems which are obtained when interactions between individual are included. The inclusion of interactions between individuals within the PTW model is actually work in progress.

From the biological viewpoint, one should not restrict the content of the paper to the sole expression of the diffusion coefficient. Indeed, once interactions between individuals will be included in the PTW model, it is not clear at all that the explicit computations which led to this expression will still be tractable. In the absence of an explicit solution of the stochastic differential system, there is little grasp to get information about the large scale behaviour of the system. By contrast, the diffusion approximation approach gives a systematic tool to study the large scale behavior of such systems, in all kinds of situations, be they linear or nonlinear. By its flexibility and its versatility, the diffusion approximation approach is the method of choice to study these problems.

One of the most popular models to describe fish behavior is the discrete Couzin-Vicsek algorithm (CVA) [START_REF] Aldana | Phase transitions in self-driven many-particle systems and related non-equilibrium models: a network approach[END_REF][START_REF] Couzin | Collective Memory and Spatial Sorting in Animal Groups[END_REF][START_REF] Grégoire | Onset of collective and cohesive motion[END_REF][START_REF] Vicsek | Novel type of phase transition in a system of self-driven particles[END_REF] (see also [START_REF] Aoki | A simulation study on the schooling mechanism in fish[END_REF][START_REF] Brillinger | Employing stochastic differential equations to model wildlife motion[END_REF][START_REF] D'orsogna | Self-propelled particles with soft-core interactions: patterns, stability and collapse[END_REF][START_REF] Mogilner | Mutual interactions, potentials, and individual distance in a social aggregation[END_REF][START_REF] Parrish | Traffic rules of fish schools: a review of agent-based approaches[END_REF][START_REF] Parrish | Self-organized fish schools: an examination of emergent properties[END_REF] for related models). For a large scale modeling of fish behavior, it is efficient to look at continuum models, which use macroscopic variables such as mean density, mean velocity and so on. Several such models based on phenomenological observations, exist (see e.g. [START_REF] Edelstein-Keshet | Mathematical models of swarming and social aggregation[END_REF][START_REF] Mogilner | A non-local model for a swarm[END_REF][START_REF] Topaz | Swarming patterns in a two-dimensional kinematic model for biological groups[END_REF][START_REF] Topaz | A nonlocal continuum model for biological aggregation[END_REF]). Several attempts to derive continuum models from the CVA model are also reported in the literature [START_REF] Kulinskii | Hydrodynamic model for a system of self-propelling particles with conservative kinematic constraints[END_REF][START_REF] Ratushnaya | Collective behaviour of self propelling particles with kinematic constraints ; the relations between the discrete and the continuous description[END_REF][START_REF] Ratushnaya | Hydrodynamic model for the system of self propelling particles with conservative kinematic constraints; two dimensional stationary solutions[END_REF]. In [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF][START_REF] Degond | Macroscopic limit of self-driven particles with orientation interaction[END_REF], a derivation of a continuum model from a kinetic version of the CVA model is proposed. However, few Individual Based Models for fish have been validated against experimental data with a comparable care as in [START_REF] Gautrais | Analyzing fish movement as a persistent turning walker[END_REF] for the PTW process. As such, the continuum model derived in this paper has a firm experimental basis, although further work needs certainly to be done to fully validate its biological foundations. Additional references on swarm aggregation and fish schooling can be found in [START_REF] Camazine | Self-Organization in Biological Systems[END_REF]. Among other types of animal societies, ants have been the subject of numerous studies and the reader can refer (among others) to [START_REF] Jost | From individual to collective ant displacements in heterogenous environments[END_REF][START_REF] Theraulaz | Spatial patterns in ant colonies[END_REF], and references therein.

The derivation of macroscopic models from particle behavior has been initiated by the seminal works of Boltzmann, and later Hilbert, Chapman and Enskog. We refer to [START_REF] Cercignani | The mathematical theory of dilute gases[END_REF] for a mathematical perspective and to [START_REF] Degond | Macroscopic limits of the Boltzmann equation: a review[END_REF] for an introduction to the subject from a modeling perspective. More recently, the derivation of macroscopic models from microscopic behavior has been very productive in other context like traffic [START_REF] Aw | Derivation of continuum traffic flow models from microscopic follow-the-leader models[END_REF][START_REF] Helbing | Traffic and related self-driven many-particle systems[END_REF] or supply-chains [START_REF] Armbruster | A model for the dynamics of large queuing networks and supply chains[END_REF]. Diffusion approximation problems for kinetic equations have been widely studied in the literature, in the context of neutron transport (see e.g. [START_REF] Bardos | Diffusion approximation and computation of the critical size[END_REF][START_REF] Bensoussan | Boundary layers and homogenization of transport processes[END_REF]), semiconductors [START_REF] Ben Abdallah | Electron transport in semiconductor superlattices[END_REF][START_REF] Degond | Diffusion approximation of a scattering matrix model of a semiconductor superlattice[END_REF][START_REF] Golse | Limite fluide des équations de Boltzmann des semiconducteurs pour une statistique de Fermi-Dirac[END_REF][START_REF] Poupaud | Diffusion approximation of the linear semiconductor equation: analysis of boundary layers[END_REF], plasmas [START_REF] Degond | Diffusion dynamics of an electron gas confined between two plates[END_REF][START_REF] Degond | Diffusion driven by collisions with the boundary[END_REF][START_REF] Degond | Existence of Solutions and Diffusion Approximation for a Model Fokker-Planck Equation[END_REF] or polymeric fluids [START_REF] Degond | Viscoelastic fluid models derived from kinetic equations for polymers[END_REF].

The outline of the paper is as follows : in section 2, the PTW process is introduced and the main results are stated. In section 3 the diffusion coefficient is obtained by direct inspection of the trajectories of the stochastic differential system. In section 4, the diffusion approximation of the forward Kolmogorov equation of the stochastic process is performed. Section 5 is devoted to proving that the trajectory method and the diffusion approximation method give rise to the same value of the diffusion coefficient. In section 6, the theoretical results are illustrated by and complemented with some numerical simulations. A conclusion is drawn in section 7. Several proofs of auxiliary results, which are inessential for the main discussion are collected in three appendices (A, B and C).

The Persistent Turning Walker model: presentation and main results

The starting point of the present work is a new model of fish motion based on experimental data taken from experiments run in La Réunion islands during years 2001 and 2002 [START_REF] Gautrais | Analyzing fish movement as a persistent turning walker[END_REF].

The studied species is a pelagic fish named Kuhlia Mugil. Its typical size ranges between 20 and 25 cm. The first experiments have been made with a single fish in a basin of 4 meters diameter during two minutes. A video records the positions of the fish every 12-th of a second (see figure 1). Then the data have been statistically analyzed and a model has been extracted [START_REF] Gautrais | Analyzing fish movement as a persistent turning walker[END_REF]. The conclusion of the statistical analysis is that the trajectories are well described by the following systems of stochastic differential equations:

d x dt = c τ (θ), (2.1
)

dθ dt = cκ, (2.2
)

dκ = -aκ dt + b dB t , (2.3)
where x = (x 1 , x 2 ) ∈ R 2 is the (two-dimensional) position vector of the (centroid of the) fish, τ (θ) = (cos θ , sin θ) is the director of the velocity vector with the angle θ ∈ Π = R/2πZ measured from the x 1 direction, κ ∈ R is the curvature of the trajectory and dB t is the standard Brownian motion. The magnitude of the velocity is constant and denoted by c > 0. The constant a is a relaxation frequency and b quantifies the intensity of the random curvature jumps. b has the dimension of 1/(L √ T ) where L and T stand for the dimensions of length and time.

The κ-dynamics is a standard Ornstein-Uhlenbeck process. The term "b dB t " models a diffusion process in curvature space while the term "-a K dt" expresses the tendency of the individual to return to a straight line trajectory. The curvature cannot increase endlessly as a consequence of the diffusion process, but rather, must relax to zero and the relaxation is stronger as the curvature gets larger. This model has been called the Persistent Turning Walker model (PTW) because it allows large excursions of the curvature towards positive or negative values, during which the spinning of the trajectory persists for a certain time.

We stress the difference with more standard diffusion processes (such as those suffered by photons in a diffusive medium), in which the Brownian motion acts on the velocity itself (or, in the case of a velocity of constant magnitude, on the angle θ). In this case, the diffusion process acts on the second derivative of the particle positions, and the associated kinetic equation is of Fokker-Planck type. This model of photon diffusion is also relevant for a certain number of animal species [START_REF] Othmer | The Diffusion Limit of Transport Equations II: Chemotaxis Equations[END_REF].

In the PTW model, the diffusion process acts on the curvature, i.e. on the third derivative of the position vector. An intuitive justification of the relevance of this model for animal behaviour is by considering the non-differentiability of the Brownian motion. Because of this feature, the photon diffusion process involves infinite second derivatives of the position, i.e. infinite forces. However, an animal body can only exert finite forces and the muscles act only in such a way that the velocity angle undergoes smooth variations. The PTW model precisely presents this feature of having smooth second order derivatives, i.e. smooth forces.

Our goal in the present work is to study the large-scale dynamics of the stochastic differential system (2.1)- (2.3). This is best done in scaled variables, where the dimensionless parameters of the model are highlighted. We use t 0 = a -1 as time unit, x 0 = ca -1 as space unit, and κ 0 = x -1 0 as curvature unit, and we introduce the dimensionless time, space and curvature as t ′ = t/t 0 , x ′ = x/x 0 and κ ′ = κ/κ 0 . For simplicity, we omit the primes. In scaled variables, the PTW model is written:

d x dt = τ (θ), (2.4
)

dθ dt = κ, (2.5
)

dκ = -κ dt + √ 2αdB t , (2.6)
where the only dimensionless parameter left is α such that

α 2 = b 2 c 2 2a 3 , (2.7)
The meaning of α 2 is the following: b/ √ a is the amplitude of a curvature change during a relaxation time a -1 , while c/a is obviously the distance travelled by the particle during this time. The product of these two quantities is dimensionless and is equal to √ 2α. It quantifies the strength of the curvature jumps relative to the other phenomena.

The individual dynamics can be translated in terms of a probability distribution f (t, x, θ, κ) d x dθ dκ of finding particles at times t with position in small neighborhoods d x dθ dκ of position x, velocity angle θ and curvature κ. The link between the individual dynamics and the evolution of the probability distribution f is given by the forward Kolmogorov equation :

∂ t f + τ • ∇ x f + κ∂ θ f -∂ κ (κf ) -α 2 ∂ κ 2 f = 0. (2.8)
This equation is an exact transcription of the individual dynamics, where the initial value f 0 at time t = 0 is given by the probability distribution of the initial conditions of the stochastic differential system (2.4)-(2.6). For more detailed considerations about the forward Kolmogorov equation and its link with stochastic differential systems, we refer the reader to [START_REF]Stochastic differential equations[END_REF][START_REF] Bass | Diffusions and elliptic operators[END_REF].

In order to capture the macroscopic dynamics, two possible routes can be taken, using either the stochastic differential system (2.4)- (2.6) or the partial differential equation (2.8). In this work, we follow both routes and verify that they lead to the same largescale behaviour. The advantage of working directly on the stochastic system is that it is simpler and it leads to explicit formulas. However, as soon as the system gets more complicated, and in particular nonlinear, explicit solutions can no longer be found and this methodology can hardly be pursued. On the other hand, the PDE approach, which, in the present case is more complicated, is also more systematic and more general. In particular, it is generally usable in the more complex nonlinear cases (see e.g. [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF][START_REF] Degond | Macroscopic limit of self-driven particles with orientation interaction[END_REF]). A particular important complex situation is the case of many interacting fish. In future work, we plan to extend the PTW model to populations of interacting fish and to use the PDE approach to extract the large-scale dynamics of the system.

From the analysis of the individual trajectories, explicit exact expressions for κ and θ in terms of stochastic integrals can be found. Unfortunately, there is no such explicit result for the position x(t), but we can calculate the first two moments of the probability distribution of x(t) explicitly, using the expressions of κ and θ. We show that the mean of the position vector stays at the origin: E{ x(t)} = (0, 0) (where E denotes the expectation over all sources of randomness, in the initial data and in the stochastic process) and that the variance grows asymptotically linearly in time. More exactly, we prove: Theorem 2.1 Under assumptions on the initial conditions that will be specified later on (see (3.1)-(3.4)), the solution of system (2.4)-(2.6) satisfies:

Var{ x(t)} t→+∞ ∼ 2D t, with D = ∞ 0 exp -α 2 (-1 + s + e -s ) ds.
(2.9)

The notation Var is for the variance over all sources of randomness. The asymptotic linear growth of the variance (2.9) suggests that the dynamics of the system is of diffusive type at large times with diffusion coefficient D. We can find an expression of D in terms of special functions. Indeed, we have Proposition 2.2 The following expression holds true:

D = e α 2 α 2 γ(α 2 , α 2 ), (2.10)
where γ(z, u) is the incomplete gamma function:

γ(z, u) = u 0 e -t t z-1 dt. (2.11)
D has the following series representation:

D = e α 2 ∞ n=0 (-1) n α 2n n! (n + α 2 ) . ( 2

.12)

To investigate the large scale dynamics of the solution of the kinetic equation (2.8) (the existence of which can be easily proved, see proposition 4.2), we need to rescale the variables to the macroscopic scale. Indeed, in eq. (2.8), all the coefficients are supposed to be of order unity. This means that the time and space scales of the experiment are of the same order as the typical time and length scales involved in the dynamics, such as, the relaxation time or the inverse of the typical random curvature excursions. Of course, in most experiments, this is not true, since the duration of the experiment and the size of the experimental region are large compared with the time and length scales involved in the dynamics.

To translate this observation, we change the space unit x 0 to a new space space unit x ′ 0 = x 0 /ε, where ε ≪ 1 is a small parameter. This induces a change of variables x ′ = εx.

We make a similar operation on the time unit t ′ 0 = t 0 /η, t ′ = ηt with η ≪ 1. Now, the question of linking η to ε is a subtle one and is largely determined by the nature of the asymptotic regime which is achieved by the system. In the present case, we expect that the asymptotic regime will be of diffusive nature, in view of theorem 2.1 and so, we will investigate the so-called 'diffusion approximation' which involves a quadratic relationship between η and ε: η = ε 2 .

For this reason, we introduce the diffusive rescaling:

t ′ = ε 2 t ; x ′ = ε x, (2.13)
and we make the following change of variable in the distribution f :

f ε (t ′ , x ′ , θ, κ) = 1 ε 2 f t ′ ε 2 , x ′ ε , θ, κ .
The scaling of the magnitude of the distribution function is unnecessary, since the problem is linear. However, it is chosen in order to preserve the total number of particles. Introducing (2.13) into (2.8) leads to the following problem for f ε :

ε∂ t f ε + τ • ∇ x f ε + 1 ε [ κ∂ θ f ε -∂ κ (κf ε ) -α 2 ∂ κ 2 f ε ] = 0 (2.14)
In order to analyze the large-scale dynamics of (2.14), we need to investigate the limit ε → 0. We show that f ε converges to an equilibrium distribution function (i.e. a function which cancels the O(ε -1 ) term of (2.14)) f 0 which depends parametrically on the particle density n 0 (x, t) and n 0 evolves according to a diffusion equation. More precisely, we prove:

Theorem 2.3 Under hypothesis 4.1 on the initial data to be precised below, the solution f ε of (2.14) converge weakly in a Banach space also to be specified below, (see (4.21)) X:

f ε ε→0 ⇀ n 0 M(κ) 2π in X weak star, (2.15)
where M is a Gaussian distribution of the curvature with zero mean and variance α 2 (see 4.4) and n 0 = n 0 (x, t) is the solution of the system:

∂ t n 0 + ∇ x • J 0 = 0, (2.16) J 0 = -D∇ x n 0 , (2.17)
where the initial datum n 0 0 and the diffusion tensor D will be defined later on (see (4.31) and (4.28) respectively).

The following theorem connects the two methods by showing that the tensor D is related to D given by (2.9): Theorem 2.4 The tensor D defined by (4.28) satisfies :

D = D 2 Id, (2.18)
where D is given by (2.9) and Id denotes the 2 × 2 identity tensor.

This theorem confirms that the trajectory method and the asymptotic PDE method are equivalent. The factor 2 between the two coefficients comes from the dimension of the problem. Indeed, D is the variance of

|x| 2 = |x 1 | 2 + |x 2 | 2 while D is the diffusion coefficient in a given direction.

Large-scale dynamics of the PTW model by the trajectory method

In this section, we want to show theorem 2.1 and proposition 2.2. We first specify the initial conditions. First, we fix the starting point of the particle at the origin :

x(t = 0) = (0, 0). (3.1)
We suppose that the initial velocity angle is uniformly distributed on the one-dimensional sphere, i.e. :

dP {θ| t=0 = θ} = dθ 2π . (3.2)
For the curvature, we make the following observation: eq. (2.6) predicts that the process κ(t) converges exponentially fast to its stationary state, which is a Gaussian distribution with zero mean and variance equal to α 2 [START_REF]Stochastic differential equations[END_REF]. We denote such a Gaussian distribution by N (0, α 2 ). For this reason, we suppose:

dP {κ| t=0 = κ} = N (0, α 2 )(κ). (3.
3)

The last hypothesis on the initial conditions is the following:

The processes θ(t = 0), κ(t = 0) and B t are independents. (

We stress that this choice of initial conditions is for simplicity only. Completely arbitrary initial conditions would lead to the same large time behaviour, but the computation would be slightly more complicated. Since we are mainly interested in the explicit expression of D, a choice of initial conditions which simplifies the calculations is legitimate.

We begin by proving the following proposition:

Proposition 3.1 The solution of the stochastic differential equation (2.4)-(2.6) with initial condition given by (3.1)-(3.4) satisfies :

E{ x(t)} = (0, 0) , ∀ t ≥ 0, (3.5) Var{ x(t)} = 2 t s=0 (t -s) exp -α 2 -1 + s + e -s ds. (3.6)
To prove this proposition, we first establish explicit formulae for the solutions of (2.5) and (2.6). The proof is deferred to appendix A.

Lemma 3.2

The solution of the stochastic differential system (2.5), (2.6) with initial conditions (3.2)-(3.4) is given by:

θ(t) = θ 0 + κ 0 -κ(t) + √ 2αB t , (3.7)
κ(t) = e -t κ 0 + √ 2αe -t t 0 e s dB s . (3.8)
Additionally,

θ(t) = θ 0 + K t 0 , (3.9)
where K t 0 is a Gaussian random variable independent of θ 0 with zero mean and variance β 2 t given by :

β 2 t = Var{K t 0 } = 2α 2 (-1 + t + e -t ). (3.10)
Proof of proposition 3.1: Using Lemma 3.2, we can compute the first two moments of x(t). Let us start with the computation of the mean. If we write x(t) = (x 1 (t) , x 2 (t)), we have :

x 1 (t) = t 0 cos θ(s) ds , x 2 (t) = t 0 sin θ(s) ds,
and, computing the mean :

E{x 1 (t)} = E t 0 cos θ(s) ds = t 0 E {cos θ(s)} ds.
Now, we can develop θ(s) using (3.9):

E {cos θ(s)} = E {cos(θ 0 + K s 0 )} = E {cos θ 0 cos K s 0 -sin θ 0 sin K s 0 } .
By the independence of θ 0 and K s 0 we finally have:

E {cos θ(s)} = E{cos θ 0 }E{cos K s 0 } -E{sin(θ 0 )}E{sin K s 0 } = 0,
since the expectations of cos θ 0 and sin θ 0 over the uniform probability distribution on θ 0 are zero. Finally, we have E{x 1 (t)} = 0, and similarly for x 2 . This proves (3.5). Now for the variance of x(t), we write:

Var{ x(t)} = E{x 2 1 (t) + x 2 2 (t)} = 2E{x 2 1 (t)}. (3.11)
by the isotropy of the problem. Then,

E{x 2 1 (t)} = E t 0 cos θ(s) ds 2 = t 0 t 0 E{cos θ(s) cos θ(u)} dsdu = 2 t 0 du u 0 ds E{cos θ(s) cos θ(u)}.
Since u ≥ s, we can write θ(u) as follows :

θ(u) = θ 0 + s 0 κ(z) dz + u s κ(z) dz = θ 0 + K s 0 + K u s
where K s 0 and K u s are Gaussian random variables independent of θ 0 with zero mean and variances β 2 s and β 2 u-s respectively, thanks to (3.10). Then, using standard identities for trigonometric functions, we get

E{cos θ(s) cos θ(u)} = E{cos(θ 0 + K s 0 ) cos(θ 0 + K s 0 + K u s )} = 1 2 (cos(2θ 0 + 2K s 0 + K u s ) + cos(-K u s )).
But since θ 0 is independent of K s 0 and K u s we have E{cos(2θ 0 + 2K s 0 + K u s )} = 0 since the mean of a cos(θ 0 + C) over the uniform distribution of θ 0 is zero whatever the value of C. Then :

E{cos(θ(s)) cos(θ(u))} = 1 2 E{cos (-K u s )} = 1 2 R cos(y) 1 √ 2πβ u-s e -y 2 2β 2 u-s dy = 1 2 e -1 2 β 2 u-s .
Indeed, an elementary computation shows that for any Gaussian random variable Z with zero mean and variance σ 2 , one has

E{cos(Z)} = exp(-σ 2 /2).) (3.12)
Thus,

E{x 1 (t) 2 } = t u=0 u s=0 exp -α 2 -1 + |u -s| + e -|u-s| dsdu.
Using the change of unknowns w = us and y = u and inverting the order of integration we find :

E{x 1 (t) 2 } = t w=0 (t -w) exp -α 2 -1 + w + e -w dw.
Using (3.11), we finally find (3.6), which ends the proof of the proposition.

In order to prove 2.1, we investigate the behavior of the variance Var{ x(t)} (given by (3.6)) when t → +∞.

End of proof of Theorem 2.1: we write, thanks to (3.6):

Var{ x(t)} -2Dt = -2 t s=0 se -α 2 (-1+s+e -s ) ds -2 ∞ s=t te -α 2 (-1+s+e s ) ds.
We have to show that the difference is bounded independently of t. For the first term, we have:

t s=0 se -α 2 (-1+s+e -s ) ds ≤ t 0 e α 2 se -α 2 s ds,
and integrating by parts, we find :

t s=0 se -α 2 (-1+s+e s ) ds ≤ e α 2 α 2 -te -α 2 t - e -α 2 t α 2 + 1 α 2 ≤ C 1 .
For the second term, we have:

∞ s=t te -α 2 (-1+s+e s ) ds ≤ t ∞ t e α 2 e -α 2 s ds ≤ te α 2 e -tα 2 α 2 ≤ C 2 .
This proves that the difference is Var{ x(t)} -2Dt is bounded independently of t and completes the proof.

We now prove Proposition 2.2 which gives an explicit approximation of the diffusion coefficient. This approximation is useful for practical simulations.

Proof of Proposition 2.2: the change of variables t = α 2 e -s in the integral (2.9) leads to (2.10). The series representation (2.12) follows from a similar series representation of the incomplete gamma function (see e.g. formula (8.354) of [START_REF] Gradshteyn | Tables of integrals, series and products[END_REF]). The series representation can also be found by expanding the exponential in the integral (2.11) in power series. This point is left to the reader.

4 Large-scale dynamics of the PTW model through the diffusion approximation of the associated kinetic equation

Formal asymptotics

In this section, for the reader's convenience, we give a formal proof of theorem 2.3. We write (2.14) as follows:

ε∂ t f ε + τ • ∇ x f ε + 1 ε Af ε = 0 (4.1)
where we define the operator A acting on functions u(θ, κ) as follows:

Au = κ∂ θ u -∂ κ (κu) -α 2 ∂ κ 2 u, (4.2)
The formal investigation of the limit ε → 0 usually starts by considering the Hilbert expansion (see e.g. [START_REF] Degond | Macroscopic limits of the Boltzmann equation: a review[END_REF] for the general theory or [START_REF] Degond | Existence of Solutions and Diffusion Approximation for a Model Fokker-Planck Equation[END_REF] for an application in the context of Fokker-Planck equations):

f ε = f 0 + εf 1 + O(ε 2 ), (4.3)
with f k being independent of ε and inserting it into (4.1). Then, collecting all the the terms of comparable orders with respect to ε, we are led to a sequence of equations. The first one, corresponding to the leading O(ε -1 ) term is Af 0 = 0, which means that f 0 lies in the kernel of A. In section 4.3, we show that the kernel of A is composed of functions of the form f 0 (t, x, θ, κ) = n 0 (t, x)M(κ)/(2π) where M(κ) is a normalized Gaussian with zero mean and variance α 2 :

M(κ) = 1 √ 2πα 2 e -κ 2 2α 2 , (4.4)
and n 0 (t, x) is a function still to be determined.

In order to determine n 0 , we first integrate (4.1) with respect to (θ, κ) ∈ Π × R and use that Au dθ dκ = 0. Defining the density n ε (t, x) and the flux J ε (t, x) by

n ε (t, x) = θ,κ f ε dκ dθ, J ε (t, x) = θ,κ f ε ε τ (θ) dκdθ. (4.5)
we find:

∂ t n ε + ∇ x • J ε = 0. (4.6)
We note that this continuity equation is valid for all values of ε. Then, letting ε → 0, we formally have n ε → n 0 . If we prove that J 0 given by (2.17) is the limit of J ε , as ε → 0, then, we can pass to the limit in (4.6) and find (2.16). System (2.16) and (2.17) is a diffusion system, which completely determines n 0 (t, x), given its initial datum n 0 0 ( x). Here, for simplicity, we assume that the initial datum for (4.1) is of the form f ε (0, x, θ, κ) = n 0 ( x)M(κ)/(2π) and the resulting initial condition for n 0 is therefore n 0 0 = n 0 (in this formal convergence proof, we admit the functions and the convergences are as smooth as required).

So, the only points left in the proof are the existence of a limit for J ε and the validity of (2.17) for J 0 . Note that the existence of a limit is not obvious because of the factor ε at the denominator of the integral (4.5) defining J ε . To prove that the limit exists, we use the Hilbert expansion (4.3) again and compute f 1 . Since A is linear, collecting the terms of order O(ε 0 ) leads to:

-Af 1 = M(κ) 2π τ • ∇ x n 0 . (4.7)
Again using the linearity of A and the fact that it operates only with respect to the (θ, κ) variables, we can write the solution of (4.7) as

f 1 = -χ • ∇ x n 0 , where χ = (χ 1 , χ 2 ) is a solution of the problem A χ = M(κ) 2π τ . (4.8)
This equation must be understood componentwise (i.e χ 1 is associated with τ 1 = cos θ and χ 2 with τ 2 = sin θ). Since the right-hand side of (4.8) has zero average with respect to (θ, κ), proposition 4.5 below shows that it has a unique solution, up to an element of the kernel of A. We can single out a unique solution by requesting that χ has zero average with respect to (θ, κ) as well. Then, all solutions f 1 to (4.7) can be written as

f 1 = -χ • ∇ x n 0 + n 1 (t, x)M(κ)/(2π) , (4.9)
where the second term of (4.9) is an arbitrary element of the kernel of A. We shall see that the determination of n 1 is unnecessary. Now, inserting the Hilbert expansion (4.3) into the integral (4.5) defining J ε , we find:

J ε (t, x) = 1 ε θ,κ f 0 τ (θ) dκdθ + θ,κ f 1 τ (θ) dκdθ + O(ε) = 0 + θ,κ f 1 τ (θ) dκdθ + O(ε), (4.10)
because f 0 is independent of θ and τ (θ) dθ = 0. Therefore, J ε has a limit when ε → 0 and this limit is given by

J 0 (t, x) = θ,κ f 1 τ (θ) dκdθ, (4.11)
To compute J 0 we insert expression (4.9) into (4.11) and find

J 0 (t, x) = θ,κ (-χ • ∇ x n 0 + n 1 M/(2π)) τ (θ) dκdθ, (4.12)
The second term vanishes and the first one can be written

J 0 (t, x) = - θ,κ τ ⊗ χ dθ dκ ∇ x n 0 , (4.13)
which is nothing but formula (2.17) with the diffusivity tensor D given by (4.28). This shows the formal convergence of the solution of the Fokker-Planck equation (2.14) to that of the diffusion system (2.16), (2.17). Now, to make this proof rigorous, we need to justify all the formal convergences. In the framework of the Hilbert expansion, this requires to work out the regularity of the various terms of the expansion. This is doable and actually leads to stronger convergences than the one we are going to prove, but this is a bit technical (see e.g. [START_REF] Degond | Existence of Solutions and Diffusion Approximation for a Model Fokker-Planck Equation[END_REF]).

What we are going to do instead is proving a convergence result in a weaker topology without using the Hilbert expansion technique. The method is close to the so-called moment method, which consists in integrating the equation against suitable test functions. This convergence proof is developed in section 4.3, but before that, we state an existence result for the original Fokker-Planck equation (4.1).

Functional setting and existence result

We define the differential operator D acting on smooth functions f (κ) by :

Df = ∂ κ (κf ) + α 2 ∂ κ 2 f. (4.14)
We state some properties of D, the proofs of which are easy and left to the reader. We recall that M(κ) denotes the normalized Gaussian with zero mean and variance α 2 (4.4).

Proposition 4.1 Let f and g be smooth functions decreasing at infinity. The following identities hold true:

Df = α 2 ∂ ∂κ M ∂ ∂κ f M , (4.15) R Df g dκ M = -α 2 R M ∂ κ f M ∂ κ g M dκ = R f Dg dκ M , (4.16)
R Df f dκ M = -α 2 R M ∂ κ f M 2 dκ ≤ 0, (4.17
)

Df = 0 ⇔ ∃c ∈ R, f = cM. (4.18)
The first identity translates the fact that M is the stationary measure of the Ornstein-Uhlenbeck process. The second one that D is formally self-adjoint with respect to the measure dκ/M. The third one shows that D is dissipative. The same inequality holds with any non-decreasing function η(f ), indeed, D(f ) η(f ) M -1 dκ ≤ 0. If η is the logarithm function, the corresponding quantity would be the relative entropy dissipation of f with respect to M. Entropy plays an important role in kinetic theory (see [START_REF] Cercignani | The mathematical theory of dilute gases[END_REF] for a review). Finally, the last quantity states that the kernel of D is one-dimensional and spanned by M.

Proposition 4.1 shows that the natural L 2 norm associated with this operator has a weight M -1 and that the natural H 1 semi-norm is given by the right-hand side of (4.17). This motivates the introduction of the following functional spaces, endowed with their naturally associated Hilbert structures and norms:

H = {u : Π θ × R κ → R / θ,κ |u(θ, κ)| 2 dθdκ M < +∞}, (4.19)
V = u ∈ H / Π,R M ∂ κ u M 2 dκdθ < +∞ , (4.20)
L 2 M = L 2 (R 2 x , H), X = L 2 ([0, T ] × R 2 x , V ). ( 4.21)
Identifying H with its dual, with have a Hilbertian triple V ⊂ H ⊂ V ′ , where V ′ is the dual of V and all injections are continuous. They are not compact because V does not bring any regularity with respect to θ.

The existence proof follows closely the existence proof of [START_REF] Degond | Global Existence of Solutions for the Vlasov-Fokker-Planck Equation in 1 and 2 Space Dimensions[END_REF] (see appendix A of this reference) and for this reason, is omitted (see also [START_REF] Degond | Existence of Solutions and Diffusion Approximation for a Model Fokker-Planck Equation[END_REF]). The proof relies on an existence theorem due to J. L. Lions [START_REF] Lions | Equations différentielles opérationnelles et problèmes aux limites[END_REF]. Proposition 4.2 Let ε > 0. We assume that f 0 belongs to L 2 M defined by (4.21). Then there exists a unique solution f ε to (2.14) with initial datum f 0 in the class of functions Y defined by :

Y = f ∈ X / ∂ t f + ε -1 τ • ∇ x f + ε -2 κ∂ θ f ∈ X ′ .
Moreover, we have the inequality for any T > 0:

||f ε (T )|| 2 L 2 M + α 2 ε 2 T 0 x,θ,κ M ∂ κ f ε M 2 dκ dθ d x dt = ||f ε (0)|| 2 L 2 M . (4.22)
Estimate (4.22) is obtained via a Green and a trace formula for functions belonging to Y which can be deduced from the one proved in [START_REF] Degond | Global Existence of Solutions for the Vlasov-Fokker-Planck Equation in 1 and 2 Space Dimensions[END_REF].

Rigorous asymptotics

We first study operator A given by (4.2), i.e. Af = κ∂ θ f -Df and state some properties which will be proved in appendix B. We view A as an unbounded operator on the Hilbert space H with domain D(A) given by:

D(A) = {u(θ, κ) ∈ V / Au ∈ H} .

Lemma 4.3

Operator A is maximal monotone. Moreover its kernel (or Null-space) is given by:

Ker(A) = {cM , c ∈ R}, (4.23)
with M defined by (4.4). The same lemma applies to the equation A * u = g.

As an application of this lemma, let χ be the solution of :

A χ = τ (θ) M 2π , (4.27)
with τ (θ) = (cos θ , sin θ). Since τ has zero average over θ and κ, χ is well-defined and unique thanks to Proposition 4.5. Then, we define the tensor D by:

D = θ,κ τ (θ) ⊗ χ dθdκ. (4.28)
Note that, since θ,κ τ (θ)M(κ) dθdκ = 0, it would not change the value of D to add any element of Ker(A) to χ.

Lemma 4.6 Let R denote the reflection operator u(θ, κ) → Ru(θ, κ) = u(θ, -κ). Then, χ * = R χ is the unique solution (satisfying (4.26)) of

A * χ * = τ (θ) M 2π , (4.29)
and we have

D = θ,κ τ (θ) ⊗ χ * dθdκ. (4.30)
Proof: Obviously, D commutes with R: DR = RD while κ∂ θ anticommutes with R: κ∂ θ (Ru) = -R(κ∂ θ u). Therefore, RA = A * R. Since the right-hand side of (4.27) is invariant by R, applying R to both sides of (4.27) leads to (4.29). Then, the change of variables κ ′ = -κ in the integral at the right-hand side of (4.30) shows that it is equal to D.

To study the limit ε → 0, we make the following hypothesis on the initial conditions.

Hypothesis 4.1 We suppose that the initial condition f ε 0 is uniformly bounded in L 2

M

and converges weakly in L 2 M to f 0 0 as ε → 0.

We can now prove theorem 2.3. The initial datum for the diffusion system (2.16, (2.17) will be shown to be:

n 0 t=0 = n 0 0 = θ,κ f 0 0 ( x, θ, κ) dκdθ. (4.31)
Proof of theorem 2.3: By hypothesis 4.1 inequality (4.22) implies :

||f ε (T )|| 2 L 2 M + α 2 ε 2 T 0 x,θ,κ M ∂ κ f ε M 2 dκ dθ d x dt ≤ C , (4.32) with C independent of ε. So (f ε ) ε is a bounded sequence in L ∞ (0, T, L 2 M ) and satisfies T 0 x,θ,κ M ∂ κ f ε M 2 dκ dθ d x dt ≤ Cε 2 , (4.33)
for any time interval T (by the diagonal process, we will eventually be able to take an increasing sequence of times T tending to infinity, so that the result will be valid on the whole interval t ∈ (0, ∞)). Therefore, there exists f 0 ∈ L ∞ (0, T, L 2 M ) and a subsequence, still denoted by f ε , satisfying :

f ε ε→0 ⇀ f 0 in L ∞ (0, T, L 2
M ) weak star . Furthermore, with (4.33), we deduce that f 0 = C( x, θ, t)M(κ). Then, letting ε → 0 in (2.14), we get that Af 0 = 0 in the distributional sense. This implies that C( x, θ, t) is independent of θ and we can write

f 0 (t, x, θ, κ) = n 0 (t, x) M(κ) 2π , (4.34)
the quantity n 0 (t, x) = f 0 (t, x, θ, κ) dθ dκ being the density associated with f 0 . Our next task is to show that n 0 satisfies the diffusion model (2.16), (2.17) with initial condition (4.31). We first note that f ε is a week solution of (2.14) with initial condition f ε 0 in the following sense: f ε satisfies:

T 0 x,θ,κ f ε (-ε∂ t ϕ -τ • ∇ x ϕ + 1 ε A * (ϕ)) dκdθd x M dt = ε x,θ,κ f ε 0 ϕ t=0 dκdθd x M . (4.35)
for all test functions ϕ in the space

C 2 c ([0, T ) × R 2 x × Π θ × R κ ) of twice continuously differentiable functions with compact support in [0, T ) × R 2 x × Π θ × R κ .
Again, the trace at t = 0 has a meaning, thanks to a trace formula for functions in Y which is proven in [START_REF] Degond | Global Existence of Solutions for the Vlasov-Fokker-Planck Equation in 1 and 2 Space Dimensions[END_REF].

We recall the definition of the flux (4.5). We prove that J ε has a weak limit as ε → 0. To this aim, in the weak formulation (4.35), we take as a test function ϕ = φ(t, x) • χ * (θ, κ) with χ * the auxiliary function defined as the solution to (4.29) and φ is a smooth compactly supported vector test function of ( x, t). Although ϕ does not have a compact support, a standard truncation argument (which is omitted here) can be used to bypass this restriction. This allows us to write:

T 0 x,θ,κ [f ε (-ε∂ t -( τ • ∇ x ))( φ • χ * ) + 1 ε f ε τ M 2π • φ] dκdθd x M dt = ε x,θ,κ f ε 0 φ t=0 • χ dκdθd x M .
Taking the limit ε → 0, we find :

lim ε→0 T 0 x J ε • φ d xdt = 2π T 0 x,θ,κ f 0 ( τ • ∇ x )( φ • χ * ) dκdθd x M dt = T 0 x n 0 ∇ x • θ,κ χ * ⊗ τ dκdθ T φ d xdt,
where the exponent T denotes the transpose of a matrix. Using (4.30) and taking the limit ε → 0 shows that J ε converges weakly (in the distributional sense) towards J 0 satisfying

T 0 x J 0 • φ d xdt = T 0 x n 0 ∇ x • (D) T φ d xdt.
This last equation is the weak form of eq. (2.17). Finally, to prove (2.16), we apply the weak formulation (4.35) to a test function of the form ϕ = φ(t, x)M(κ), where again, φ(x, t) is a scalar, smooth and compactly supported test function of ( x, t) in R 2 × [0, T ). This gives :

- T 0 x,θ,κ f ε ((ε∂ t + ( τ • ∇ x ))ϕ) dκdθd xdt = ε x,θ,κ f ε 0 φ t=0 dκdθd x.
Dividing by ε and taking the limit ε → 0, we get :

- T 0 x (n 0 ∂ t ϕ + J 0 • ∇ x φ) d xdt = x n 0 0 φ t=0 d x,
where n 0 0 is defined by (4.31). This last equation is exactly the weak formulation of equation (2.16), with initial datum n 0 0 . This concludes the proof.

Equivalence of the two methods

In this section, we show that both methods lead to the same value of the diffusion coefficient (theorem 2.4). The first step is to show that we can approximate the solution of equation (4.24) by the solution of the associated evolution equation. More precisely, in appendix C, we prove the following lemma : Lemma 5.1 Let g in H satisfying (4.25) and u ∞ in D(A) be the solution of (4.24) satisfying (4.26). Let u 0 ∈ D(A) satisfying (4.26). Then, the solution u(t) of the evolution problem:

∂ t u = -Au + g , u t=0 = u 0 , (5.1)
weakly converges to u ∞ in H as t tends to ∞.

With this lemma we can explicitly calculate the tensor D and prove the theorem 2.4 :

Proof of theorem 2.4: Let χ(t) be the solution of

∂ t χ = -A χ + τ (θ) M(κ) 2π , χ(t = 0) = 0. (5.2)
Thanks to Lemma 5.1, χ(t) weakly converges to χ in H when t → ∞. It follows tht :

κ,θ χ(t) ⊗ τ dκdθ t→+∞ -→ κ,θ χ ⊗ τ dκdθ. (5.3)
Let us consider the first component of χ(t), which we denote by u(t) and the integrals κ,θ u(t) cos θ dκdθ and κ,θ u(t) sin θ dκdθ. Because u satisfies (5.2), it admits the following representation (see the proof of Lemma 5.1):

u(t) = t 0 T s cos θ M(κ) 2π ds,
where T t is the semi-group generated by -A (see [START_REF] Pazy | Semigroups of linear operators and applications to partial differential equations[END_REF]). With this expression, we evaluate the integral of u(t) against cos θ :

κ,θ u(t) cos θ dκdθ = κ,θ t 0 T s cos θ M(κ) 2π ds cos θ dκ dθ = t 0 κ,θ cos θ M(κ) 2π T * s (cos θ) dκ dθ ds,
where T * is the adjoint operator of T in L 2 (θ, κ) generated by -A * , where

A * (f ) = -κ∂ θ f + κ∂ κ f -α 2 ∂ κ 2 f.
Note that we are referring here to the adjoint in the standard L 2 sense and not in the weighted space H. This is why A * does not coincide with A * defined in Lemma 4.4. The semi-group T * t admits a probabilistic representation: for all regular functions f (θ, κ) T * t (f )(θ, κ) = E{f (θ(t), κ(t))|θ 0 = θ, κ 0 = κ}, where (κ(t), θ(t)) is the solution of the stochastic differential equation (2.5), (2.6). Using this representation, we have :

κ M(κ) T * s (cos θ) dκ = κ M(κ) E{cos θ s |θ 0 = θ, κ 0 = κ} dκ = E{cos θ s |θ 0 = θ, κ 0 = Z},
where Z is a random variable independent of B t with density M. Using lemma 3.2, we have :

E{cos θ s |θ 0 = θ, κ 0 = Z} = E{cos(θ + Y s )},
with Y s a Gaussian random variable with zero mean and variance β 2 s given by (3.10). Then:

E{cos(θ + Y s )} = E{cos θ cos Y s -sin θ sin Y s } = cos θ E{cos Y s },
because the density of Y s is even and implies that E{sin Y s } = 0. Finally using (3.12), we have:

κ M(κ) T * s (cos θ) dκ = cos θ e -β 2 s 2 .

20

Then, the first integral is given by:

κ,θ u(t) cos θ dκ dθ = t 0 θ cos θ 2π cos θ e -β 2 s 2 dθ ds = t 0 1 2 e -β 2 s 2 ds.
We can proceed similarly to evaluate the integral of u(t) against sin(t). This gives:

κ,θ u(t) sin θ dκdθ = t 0 θ cos θ sin θ 2π E{cos Y s } dθds = 0.
It remains to evaluate the integrals involving the second component of vector χ(t) which we denote by v(t). By the same method as for u(t), we get :

κ,θ v(t) cos θ dκ dθ = 0 and κ,θ v(t) sin θ dκ dθ = t 0 1 2 e -β 2 s 2 dt.
Collecting these formulae, we can write:

κ,θ χ(t) ⊗ τ dκdθ = D(t) 2 Id, with D(t) =
t 0 e -α 2 (-1+u+e -u ) du. Taking the limit t → +∞ and using equation (5.3), shows that (2.18) holds true and completes the proof of the theorem.

Numerical simulation

We simulate individual trajectories satisfying equation (2.4)-(2.6) with initial conditions given by (3.1)- (3.4). If we fix a time step ∆t, using (3.7), (3.8), we have :

κ (n+1)∆t = γκ n∆t + G (n+1) θ (n+1)∆t = θ 0 + κ 0 -κ (n+1)∆t + √ 2αB (n+1)∆t .
with γ = e -∆t and G (n+1) a Gaussian random variable with zero mean and variance 2α 2 1e -2∆t independent of κ n∆t . With this formula, we can simulate recursively the process (κ n∆t , θ n∆t ) n exactly (in the sense that it has the same law as the exact solution).

To generate the Brownian motion, we just compute the increments B (n+1)∆t -B n∆t since they are Gaussian and independent of B n∆t . On the other hand, these increments are not independent of G (n+1) . Fortunately, we can compute the covariance matrix of the Gaussian vector (G (n+1) , B (n+1)∆t -B n∆t ) :

G (n+1) B (n+1)∆t -B n∆t ∼ N (0, C)
where N (0, C) is a two-dimensional Gaussian vector with zero mean and covariance matrix C given by:

C = 2α 2 1 -e -2∆t √ 2α(1 -e -∆t ) √ 2α(1 -e -∆t ) ∆t .
Knowing this covariance matrix, we can simulate the Gaussian vector (G (n+1) , B (n+1)∆t -B n∆t ) using the Cholesky method : we generate (X 1 , X 2 ) a vector of two independent normal law, and take √ C(X 1 , X 2 ) T as realization of the Gaussian vector. Now for the position x, since we do not have any explicit expression, we use a discrete approximation scheme of order O((∆t) 2 ). For example, the first component x 1 of x is approximated by:

x 1 ((n + 1)∆t) = x 1 (n∆t) + (n+1)∆t n∆t cos θ(s) ds ≈ x 1 (n∆t) + ∆t 2 (cos θ(n∆t) + cos θ((n + 1)∆t)).
We present four trajectories obtained with different values of the parameter α in figure 2. As the parameter α increases, the excursions towards large positive or negative curvatures become larger. As a consequence, the spinning of the trajectory around itself increases and, from almost a straight line when α = 0.1, the trajectory shrinks and looks closer and closer to a wool ball. We can figure out that the diffusion coefficient decreases as α increases which is confirmed by formula (2.9) since the function -1+s+exp(-s) ≥ 0 for s ≥ 0. Additionally, it is easily seen that D → ∞ like α -2 when α → 0 while D → 0 when α → ∞.

To illustrate theorem 2.1, we use a Monte-Carlo method to simulate the variance of the process x. We simulate N independent trajectories and we compute the variance of the sample at each time step. In figure [START_REF] Armbruster | A model for the dynamics of large queuing networks and supply chains[END_REF], we compare the result obtained with N = 2000 and the theoretical prediction given by the (3.6). The figure shows an excellent agreement between the computation and the theoretical prediction. Additionally, after an initial transient, the growth of the variance is linear, in accordance with the theoretical result (2.9). We can use the slope of the asymptotically linear part of the curve to give a numerical estimate of the diffusion coefficient D. Fro this purpose, we fit a straight line (in the mean-square sense) between times T /2 and T . We remove the data between 0 and T /2 because the initial transient is not linear and including them would deteriorate the accuracy of the measurement. We compare the slope of the fitted line with the theoretical value (2.9). We report the result of this comparison for two values of α (α = 0.1 and α = 2., with T = 1200 time units in Table 1. The approximation is quite good, with an error comprised between 2 and 3%, which can be attributed to numerical noise and to an unsufficient approximatin of the asymptotic state. To illustrate the influence of the initial transient, we take T = 120 time units in the case α = 0.1 and report the result in Table 2. There, the approximation is quite poor, because the asymptotic state has not yet been reached.

In order to illustrate theorem 2.3, we plot the spatial density n(t, x) of the distribution f (t, x, θ, κ) using a Monte-Carlo algorithm for α = 2 and T = 30 time units on figure 4. We see that the density has the Gaussian shape of the solution of a diffusion equation, in accordance with the prediction of the theorem. To make a more quantitative comparison, we compare it with the asymptotic prediction, i.e. the solution of the diffusion equation (2.16), (2.17), by computing the difference in L 1 norm. The results are reported in figure 5. We plot the L 1 norm of the difference for α = 1. and for four values of the parameter ε : 1, 1 2 1 5 and 1 10 . As expected, the agreement is better as ε is smaller. However, at large times, all solutions are eventually close to the solution of the diffusion equation. Roughly speaking, the time at which the solution of the diffusion equation starts to be a good approximation of the solution of the kinetic equation scales like ε. This means that, after an initial transient the duration of which may depend on ε, the solution is close to that of the diffusion equation, no matter the value of ε. Figure 5: L 1 norm of the difference between n ε (t, ) = f ε (t, •, θ, κ) dθ dκ and its asymptotic limit n 0 (t) as a function of time t, with α = 1 and N = 10 4 simulation particles.

Conclusion

In this paper, the large-scale dynamics of the 'Persistent Turning Walker' (PTW) model of fich behavior has been analyzed. It has been shown, by two different methods, that the large scale limit of this model is of diffusion type, and an explicit formula for the diffusion coefficient has been provided. While the direct analysis of the stochastic trajectories provides a direct route to the value of the diffusion constant, the diffusion approximation of the associated forward Kolmogorov equation, which is of Fokker-Planck type, gives a more systematic way to extend the theory to more complex nonlinear cases. Such a nonlinear situation will be encountered when, in the near future, the nonlinear interactions between the individuals will be introduced within the PTW model. We expect that, in this context, the diffusion approximation methodology will have to be exploited thoroughly to allow access to the large scale behaviour of the system.

Appendix A: proofs of section 3 Proof of Lemma 3.2: formula (3.8) is standard in the theory of Ornstein-Uhlenbeck processes [START_REF]Stochastic differential equations[END_REF]. To obtain (3.7), we integrate κ(t) with respect to time: Interchanging the order of integrations and integrating with respect to s, we deduce:

t 0 κ(s) ds = (1 -e -t )κ 0 + √ 2α t 0 (1 -e -(t-u) ) dB u .
Then we develop the integral :

t 0 κ(s) ds = (1 -e -t )κ 0 + √ 2αB t -(κ(t) -e -t κ 0 ).
This formula can be rewritten:

t 0 κ(s) ds = κ 0 -κ(t) + √ 2αB t ,
which easily leads to (3.7). We now calculate the mean and the variance of K t 0 = t 0 κ(s) ds. Since κ(s) is of zero mean, its integral K t 0 is also of zero mean: E{K t 0 } = 0. Now for the variance of K t 0 , we can write :

Var{K t 0 } = E (1 -e -t )κ 0 - √ 2αe -t t 0 e s dB s + √ 2αB t 2 .
Using that κ 0 and B s are independent, we can develop the square and get: Since (U λ ) λ is a bounded sequence in H, we can extract a subsequence (still denoted by U λ ) such that U λ ⇀ U in H weak as λ → 0. Taking the limit λ → 0 in (7.5), gives A(U) = 0. If we take the inner product of (7.5) with U λ and then pass to the limit λ → 0, we also find that U belongs to V . So Lemma 4.3 applies and gives U = cM with a constant c ∈ R. Using (4.25), we also have :

Var{K t 0 } = (1 -e -t ) 2 E{κ 2 0 } + 2α 2 E -e -t
< λU λ + AU λ , M > H =< g N λ , M > H = 0.
So κ,θ U λ dθdκ = 0 for all λ. Taking the limit λ → 0 leads to κ,θ U dθdκ = κ,θ CM(κ) dθdκ = 2πC = 0, which implies U = 0. This proves:

U λ ⇀ 0 in H weak . (7.6)
To get a contradiction, we now prove that the convergence is strong.

To this aim, we introduce a decomposition of the space H into two orthogonal subspaces. Let L be the closed subspace of H defined by :

L = {c(θ)M / c(θ) ∈ L 2 (Π θ )},
with M defined by (4.4). So H = L ⊥ ⊕ L ⊥ . We also define the orthogonal projector P of H onto L such that P f = κ f (κ, θ) dκ M. Using this projection, we decompose the sequence (U λ ) λ as follows:

U λ = c λ (θ)M + v λ , (7.7)
with v λ ∈ L ⊥ , i.e. κ v λ dκ = 0. To demonstrate that U λ λ→0 -→ 0 in H strongly, we first demonstrate that v λ λ→0 -→ 0 in H strongly. Taking the inner product of the equation satisfied by U λ (7.5) with U λ gives :

λ U λ 2 H + θ,κ M ∂ κ U λ M 2 dθdκ = 1 N λ < g, U λ > H . Since ∂ κ U λ M = ∂ κ v λ
M and U λ H = 1, we get by taking the limit λ → 0 :

θ,κ M ∂ κ v λ M 2 dθdκ λ→0 -→ 0. ( 7.8)
Now Gross inequality [START_REF] Gross | Logarithmic Sobolev inequalities and contractivity properties of semigroups[END_REF] gives, for any v ∈ V :

R ∂ κ f M 2 M dκ + R f dκ 2 ≥ R |f | 2 dκ M . (7.9)
Then, since κ v λ dκ = 0, we deduce:

R ∂ κ v λ M 2 M dκ ≥ R |v λ | 2 M dκ.
Integrating this inequality with respect to θ and using (7.8), we find:

v λ H λ→0 -→ 0, in H strong. (7.10)
To prove the convergence of c λ , we define the bounded operator T : H → L 2 (Π θ ) such that T f = κ κf dκ. Having T acting on (7.5) and taking the limit λ → 0, leads to:

T A U λ λ→0 -→ 0 in L 2 (Π θ ) strong. (7.11)
If we develop the left-hand side, we find:

T A U λ = κ κ 2 ∂ θ U λ dκ - κ κ∂ κ (κU λ ) -α 2 κ∂ κ 2 U λ dκ = κ κ 2 ∂ θ U λ dκ + κ κU λ dκ.
But using the decomposition U λ = c λ M + v λ (7.7), we have:

κ κU λ dκ L 2 (θ) = κ κv λ dκ L 2 (θ) λ→0 -→ 0.
So, (7.11) leads to :

κ κ 2 ∂ θ U λ dκ λ→0 -→ 0 in L 2 (Π θ ) strong. (7.12)
If we define h λ (θ) = κ κ 2 U λ dκ, (7.12) is equivalent to saying that ∂ θ h λ L 2 (θ) λ→0 -→ 0. Using the Poincare-Wirtinger inequality [START_REF] Brézis | Analyse fonctionnelle[END_REF], there exists a constant C 0 such that:

h λ -hλ L 2 (θ) ≤ C 0 ∂ θ h λ L 2 (θ) , (7.13)
with hλ = 1 2π 2π 0 h λ (θ) dθ. Then, we develop hλ . We get:

hλ = 1 2π 2π 0 κ κ 2 U λ dκdθ =< U λ , Mκ 2 > H λ→0 -→ 0 in R,
since U λ converges weakly to zero (see (7.6)). So, (7.13) leads to h λ λ→0 -→ 0 in L 2 (Π θ ) strong. If we develop h λ we find:

h λ (θ) = κ κ 2 (c λ (θ)M + v λ ) dκ = α 2 c λ (θ) + κ κ 2 v λ dκ.
Now, κ κ 2 v λ dκ converges to zero in L 2 (θ) strong because of (7.10) and we finally have :

c λ (θ) λ→0 -→ 0 in L 2 (θ) strong.
Using the convergence of c λ and v λ , we can now prove the strong convergence of U λ to 0 in H:

U λ 2 H = c λ M 2 H + v λ 2 H = c λ 2 L 2 (θ) + v λ 2 H λ→0 -→ 0,
which contradicts the fact that U λ has unit norm in H. This shows that there exists a bounded subsequence in the sequence u λ . In fact, since the same proof can be applied to any subsequence, this shows that the whole sequence u λ is bounded, but this is useless for our purpose. We conclude the proof of Proposition 4.5 as follows: there exists a subsequence u λ and a function u in H such that u λ ⇀ u in H weak. Taking the limit of (7.4) as λ → 0, we deduce that Au = g in the sense of distributions. However, since g ∈ H, eq. Au = g also holds in H. Moreover if we take the inner product of (7.4) with u λ and pass to the limit λ → 0, we find that u belongs to V . So u belongs to D(A), which ends the proof of the 'if' part of the statement.

Finally, to prove uniqueness, we just remark that, two solutions of (4.24) differ from an element of the kernel of A and we apply (4.23). This ends the proof.

Appendix C: proofs of section 5

Proof of Lemma 5.1: the proof borrows some ideas from [START_REF] Desvillettes | On Long Time Asymptotics of the Vlasov-Poisson-Boltzmann Equation[END_REF], but is simpler, due to the linear character of the problem. The difficulty is getting some compactness in time. Here, instead of considering time translates of the solution as in [START_REF] Desvillettes | On Long Time Asymptotics of the Vlasov-Poisson-Boltzmann Equation[END_REF], we will consider time integrals over a fixed interval length ∆t.

Since operator A is maximal monotone on H (see Lemma 4.3), operator -A generates a semi-group of contractions T t on H. Moreover the solution of (5.1) is given by: u(t) = T t (u 0 ) + t 0 T s (g) ds.

We define f (t) = u(t)u ∞ which satisfies :

∂ t f = -Af, f t=0 = f 0 , (7.1)
with f 0 = u 0u ∞ and κ f 0 (κ) dκ = 0. To prove the weak convergence of u(t) to u ∞ , we have to prove that f (t) converges to zero weakly in H.

To this aim, we make an orthogonal decomposition of f (t) as in the proof of Proposition 4.5: f (t) = c(t)M + v(t), with c(t) ∈ L 2 (Π θ ), v(t) ∈ H and κ v(t) dκ = 0. Taking the inner product of (7.1) with f , we get :

1 2 ∂ t f 2 H = - κ,θ α 2 M ∂ κ f M 2 dκdθ.
Using the decomposition of f (t) and noticing that ∂ κ f M = ∂ κ v M , this equality becomes: It remains to prove the convergence of c(t) to zero. We integrate (7.1) with respect to κ. This gives, using that κ M(κ) dκ = 1 and κ v(t) dκ = 0 : We fix a time interval ∆t and integrate (7.4) over this time interval. This leads to: Since c belongs to L ∞ ((0, ∞) t , L 2 (Π θ )) (see (7.2), we have t+∆t t c(s) ds ∈ L ∞ ((0, ∞) t , L 2 (Π θ )). So there exists a subsequence such that t+∆t t c(s) ds is weakly convergent in L 2 (Π θ ). Actually, (7.5) implies that there exists a constant function with respect to θ, depending on ∆t and denoted by L(∆t) such that To deduce the convergence of c(t), we have to control the derivative of c(t) in time. For this purpose, we rewrite :

1 2 ∂ t c(t) 2 L 2 + v(t)
∂ t c(t) =

Figure 1 :

 1 Figure 1: One trajectory of a Kuhlia Mugil fish

Lemma 4 . 4 Proposition 4 . 5

 4445 The adjoint A * of A in H is given by A * f = -κ∂ θ f -Df . It is a Maximal monotone operator with domain D(A * ) = D(A) and Ker(A * ) = Ker(A). Let g ∈ H. Then, there exists u ∈ D(A) such that Au = g, (4.24)if and only if g satisfies the following solvability condition:Moreover, the solution u is unique up to a constant times M. A unique solution can be singled out by presribing the condition

Figure 2 :

 2 Figure 2: Four trajectories simulated with different value of α: α = 0.1 (top left), α = 0.5 (top right), α = 1. (bottom left) and α = 2. (bottom right). The simulation is run during 120 time units with a time step dt = 0.05 time unit

Figure 3 :

 3 Figure 3: Variance of the process x(t): comparison between the numerical simulation (points) and the theoretical prediction (solid line) for two values of α : α = 0.1 (left) and α = 2 (right).

3 Figure 4 :

 34 Figure 4: The spatial distribution n for α = 2 at time T = 30 time units.

2α t 0 s 0 e

 0 ) ds = (1e -t )κ 0 + √ -s e u dB u ds.

  ∂ θ κ κv(t) dκ. (7.3)Now if we pre-multiply by κ before integrating with respect to κ, we obtain :∂ t κ κv(t) dκ = α 2 ∂ θ c(t) + ∂ θ κ κ 2 v(t) dκκ κv(t) dκ.(7.4)

κ

  κ(v(t + ∆t)v(t)) dκ = α 2 ∂Since v(t) converges to zero in H, we have, in the sense of distributions:

  + z) dκ dzds.Using again the convergence of v(t) to zero, we find :∆t c(t) t→+∞ ⇀ L(∆t), or defining the constant C = L(∆t)∆t , we have c(t) t→+∞ ⇀ C in L 2 (Π θ ) weak. To complete the proof, it remains to prove that C is equal to zero. Now, since eq. (7.1) is mass preserving i.e.:∂ t κ,θ f (t) dκdθ = -

  Af (t) dκdθ = 0, we have κ,θ f (t) dκdθ = κ,θ f (0) dκdθ = 0. Also :

  C = 0. This proves f (t) t→+∞ ⇀ 0 in H weak and completes the proof.

Table 1 :

 1 Diffusion coefficient: comparison of the numerical estimate obtained by fitting the numerical values with a straight line over the time interval [T /2, T ] with the theoretical prediction (2.9). T = 1200 units of time.

		D simulation D theoretical relative error
	α = 0.1	58.8	101	72 %

Table 2 :

 2 Same as Table1but with T = 120 units of time. The agreement is poor because the asymptotic state is not reached.

  where the Ito correction term is zero due to the fact that exp s is a deterministic process. We can simplify this expression again since E{B t B s } = min(t, s) and get:

								t			2
									e s dB s + B t	.
								0
	Let us consider the second term. By Ito's formula, we have
			t	2		t		2			t
	E	-e -t	0	e s dB s + B t	= e -2t E	0	e s dB s	-2e -t E B t	0	e s dB s + E B 2 t
					t
		= e -2t e E 0 -e -t t 0 e s dB s + B t 2 = 1 -e -2t 2	-2e -t e t t -	0	t	e s s ds + t
					=	1 -e -2t

2s ds -2e -t E B t (e t B t -t 0 e s B s ds + t, 2 -2 1e -t + t.

Using also that E{κ 2 0 } = α 2 , the variance of K t 0 is written:

Developing and simplifying the expression, we find (3.10), which ends the proof.

Appendix B: proofs of section 4

Proof of Lemma 4.3: Let u ∈ D(A). Then, κ∂ θ u ∈ V ′ and Lemma A1 of [START_REF] Degond | Global Existence of Solutions for the Vlasov-Fokker-Planck Equation in 1 and 2 Space Dimensions[END_REF] shows that the Green formula for functions u ∈ V such that κ∂ θ u ∈ V ′ is legitimate. Therefore, taking the inner product of A(u) against u, we find:

So, A is a monotone operator in H. To show that A is maximal monotone, we prove that for any g ∈ H, there exists u ∈ D(A) such that :

Taking the inner product of (7.2) against a test function ϕ in the space D(Π θ × R κ ) of infinitely differentiable and compactly supported functions on Π θ × R κ leads to the variational problem :

Again, the same theory as in the appendix A of [START_REF] Degond | Global Existence of Solutions for the Vlasov-Fokker-Planck Equation in 1 and 2 Space Dimensions[END_REF] (based the result by J. L. Lions in [START_REF] Lions | Equations différentielles opérationnelles et problèmes aux limites[END_REF]) applies to prove the existence of a solution to (7.3) with u in V such that κ∂ θ u ∈ V ′ . From there, it immediately follows that u ∈ D(A).

It is immediate to see that any function of the form u(θ, κ) = CM(κ) for any constant C belongs to the kernel of A. Conversely, suppose that u ∈ Ker A. Then, by (7.1), there exists a function C(θ) ∈ L 2 (Π) such that u(θ, κ) = C(θ)M(κ). But again, A(u) = 0 implies that κ ∂ θ C(θ) M = 0. So C(θ) is a constant, which proves (4.23).

Proof of Proposition 4.5: the 'only if' part of the theorem is obvious since, using Green's formula (again, obtained by adapting that of appendix B of [START_REF] Degond | Global Existence of Solutions for the Vlasov-Fokker-Planck Equation in 1 and 2 Space Dimensions[END_REF], we have Au dθ dκ = 0. To prove the 'if' part, we borrow a method from (for instance) [START_REF] Castella | Diffusion dynamics of classical systems driven by an oscillatory force[END_REF]. To find a solution to (4.24), we look at a perturbed equation :

with λ > 0. Since A is maximal monotone in H (Lemma 4.3), eq. ( 7.4) admits a solution u λ for all positive λ ( [START_REF] Brézis | Analyse fonctionnelle[END_REF]). To prove the existence of a solution to (4.24), we want to extract a subsequence, still denoted by (u λ ) which converges weakly in H. For this purpose, it is enough to show that there exists a bounded subsequence. We proceed by contradiction, supposing that the (full) sequence N λ = u λ H λ→0 → +∞. We define U λ = u λ N λ . U λ satisfies U λ H = 1 for all λ and λU λ + AU λ = g N λ .

(7.5)