Julien Barré
email: julien.barre@univ-orleans.fr

Dan Crisan
email: d.crisan@imperial.ac.uk

Thierry Goudon
email: thierry.goudon@inria.fr

Two-dimensional pseudo-gravity model

Keywords: Attractive forces. Convection-diffusion equations. Mean field regime. 2010 MSC Subject Classification. 82C70, 70F45, 35Q35

We analyze a simple macroscopic model describing the evolution of a cloud of particles confined in a magneto-optical trap. The behavior of the particles is mainly driven by self-consistent attractive forces. In contrast to the standard model of gravitational forces, the force field does not result from a potential; moreover, the non linear coupling is more singular than the coupling based on the Poisson equation. We establish the existence of solutions, under a suitable smallness condition on the total mass, or, equivalently, for a sufficiently large diffusion coefficient. When a symmetry assumption is fulfilled, the solutions satisfy strengthened estimates (exponential moments). We also investigate the convergence of the N -particles description towards the PDE system in the mean field regime.

Introduction

This work is concerned with a simple mathematical model describing anisotropic magneto-optical traps (MOT). In these devices, clouds of atoms are held together at very low temperatures through the action of well tuned lasers. These lasers induce on each atom an external space dependent confining force, as well as a friction: these effects are responsible for the trapping and cooling of the atoms. The lasers also create effective interaction forces between the atoms. The precise description of these forces involves a full description of the laser field and its coupling with the atoms. The following simplification, while probably not always quantitatively accurate, is customary since the pioneering article [START_REF] Sesko | Behavior of neutral atoms in a spontaneous force trap[END_REF]: the interaction forces are divided into i) a repulsive force due to multiple diffusion of photons, which is usually approximated by a Coulomb force (predicted in [START_REF] Sesko | Behavior of neutral atoms in a spontaneous force trap[END_REF]) and

ii) an attractive long-range force, the so-called "shadow effect" (predicted in [START_REF] Dalibard | Laser cooling of an optically thick gas: the simplest radiation pressure trap ?[END_REF]), that bears some similarity with gravity, and is the main subject of this article.

In a standard, roughly spherical, cloud, the repulsive force dominates. Nevertheless, if an external potential forces the cloud into a very elongated cigar shape, or a very thin pancake shape, the attractive force is expected to dominate, and the repulsive force may be neglected in a first approximation [START_REF] Barré | Nonequilibrium phase transition with gravitational-like interaction in a cloud of cold atoms[END_REF][START_REF] Chalony | Long-range one-dimensional gravitational-like interaction in a neutral atomic cold gas[END_REF]. This is the regime we are interested in. A typical MOT involves 10 6 to 10 10 interacting particles. Although in experiment trapping the atoms in a pancake-shaped cloud would probably contain less atoms, it is then relevant to make use of a partial differential equations describing the particles' density, instead of considering the dynamics of the individual particles. A reasonable model may be a 3D non-linear Fokker-Planck, or a McKean-Vlasov, equation. However, in order to describe the cigar-or pancake-shaped clouds observed in the experiments it makes sense to use a large scale approach, and to integrate over the small dimension(s). After some approximations, one is left with an effective 1D or 2D nonlinear partial differential equations (PDE). The 1D equation obtained this way coincides with the meanfield description a 1D damped self-gravitating system [START_REF] Chalony | Long-range one-dimensional gravitational-like interaction in a neutral atomic cold gas[END_REF] and is well-known. We thus concentrate on the 2D case. The 2D nonlinear PDE studied here has its own interest, independently of the relation with the MOT experiments: it bears some similarities with a 2D damped self-gravitating system (also known as the Smoluchowski model in astrophysics [START_REF] Chandrasekhar | Brownian motion, dynamical friction and stellar dynamics[END_REF] or the Keller-Segel chemotactic model [START_REF] Keller | Initiation of slime mold aggregation viewed as an instability[END_REF][START_REF] Keller | Model for chemotaxis[END_REF]). Therefore, a natural question is to determine whether or not singularities appear in finite time, depending on certain thresholds, as this is the case for the Keller-Segel model, see the review [START_REF] Horstmann | From 1970 until present: The Keller-Segel model in chemotaxis and its consequences[END_REF][START_REF] Horstmann | From 1970 until present: The Keller-Segel model in chemotaxis and its consequences[END_REF].

We are interested in the particle density (x, y, t) → ρ(x, y, t), which is a scalar non-negative quantity that depends on the time t ≥ 0 and space variables (x, y) ∈ R 2 . Its evolution is governed by the following non linear PDE

∂ t ρ = ∇ • D∇ρ -F [ρ]ρ , (1)
where the constant D > 0 is given and the self consistent force field

F [ρ] = F x [ρ] F y [ρ] ,
is defined by F x [ρ](x, y, t) = -sgn(x -x ′ )ρ(x ′ , y, t) dx ′ , F y [ρ](x, y, t) = -sgn(y -y ′ )ρ(x, y ′ , t) dy ′ .

(

) 2
The problem is complemented with an initial data

ρ t=0 = ρ 0 . (3)
Similar to the Keller-Segel model, the force is thus defined through a convolution formula. As a consequence of the fact that the (distributional) derivative of the function x → sgn(x) is 2δ 0 , where δ 0 is the Dirac delta distribution at 0, we observe that (mind the sign)

∇ • F [ρ] = -4ρ ≤ 0. ( 4
)
The divergence of the force field of the Keller-Segel system satisfies the same relation. However, there are crucial differences with the Keller-Segel system that make the analysis here different:

• the force does not have the potential structure ( F cannot be expressed as the gradient of a potential), and, accordingly, we cannot derive estimates related to the evolution of a potential energy,

• the convolution acts only on a single direction variable; hence we cannot expect any regularisation effect similar to the one given by the coupling of the force through the Poisson equation,

F (ε) y [ρ](x, y, t) = - sgn (ε) (y -y ′ )δ (ε) (x -x ′ )ρ(x ′ , y ′ , t) dx ′ dy ′ (6)
with

sgn (ε) (u) = 2 1 ε √ 2π u 0 e -v 2 2ε 2 dv, δ (ε) (u) = 1 2 d du sgn (ε) (u) = 1 ε √ 2π e -u 2 2ε 2 .
Denoting by ⋆ x (resp. ⋆ y ) the convolution with respect to the variable x (resp. the variable y), we observe that

F (ε) x [ρ] = T (ε) (sgn ⋆ x ρ) = sgn ⋆ x (T (ε) ρ), F (ε) y [ρ] = T (ε) (sgn ⋆ y ρ) = sgn ⋆ y (T (ε) ρ),
where T (ε) stands for the convolution with the normalized 2-d Gaussian kernel.

2. Establish a priori estimates that are uniform with respect to ε. We obtain several such estimates, typically L p and moment estimates, based on dissipative properties of the equation, at the price of assuming the diffusion coefficient D large enough. Section 2 includes these estimates.

3. Show the existence and uniqueness of solutions ρ (ε) of the regularized PDE [START_REF] Brézis | Analyse fonctionnelle. Théorie et applications[END_REF]. To this end, we employ a suitable fixed point approach, described in Section 3

4. Use the a priori estimates to prove global existence of the solution of the original equation, at least when D is large enough. We present two proofs. The first relies on quite standard compactness arguments. As mentioned above the difficulty is related to the non-linear term F [ρ]ρ and the adopted functional framework should be constructed so that the product makes sense and is stable. The second approach is more precise and establishes directly that the sequence of approximated solutions ρ (ε) ε>0 satisfies the Cauchy criterion in a certain norm. However this approach requires certain symmetry assumptions and fast enough decay of the initial state. These additional assumptions allow us to derive exponential moments, and weighted estimates on the gradient of the unknown. This analysis is detailed in Section 4.

5. Introduce a stochastic system of N ≫ 1 particles, with a regularized interaction, and prove that the empirical measure converges towards a solution of the PDE when N → ∞. Assuming that the number of particles N is proportional to e C/ε 2 with ε being the regularizing parameter, one can obtain particle approximations that are arbitrarily close to ρ, on any fixed time interval. In particular we show that one can get an upper bound for the Wasserstein distance between the particle approximation and ρ of order ε -ν , where ν > 0 is a certain constant independent of ε. Put it differently, we show that if ε is of order (log(N )) -1 2 , then the rate of convergence of the Wasserstein distance between the particle approximation and ρ is also of logarithmic order, see Theorem 5.2. The analysis is presented in Section 5. 6. Run numerical simulations using the particle representation obtained in Section 5 and compare it with the PDE method introduced in [START_REF] Carrillo | A finite-volume method for nonlinear nonlocal equations with a gradient flow structure[END_REF]. In this way we illustrate the existence results covered by Theorems 4.1 and 4.5. As we will see, the constraint on the diffusion coefficient (condition ( 12)) required for the two theorems to be valid is not optimal: the solution can apparently be global in time for other values too. We also illustrate the convergence for the particles approximation. The rate of convergence of the particle approximation as a function of N seems to be much better than that suggested by Theorem 5.2. These are covered in Section 6.

A priori estimates 2.1 Moments

Let k ∈ N, k ≥ 2. We set m k (t) = (|x| k + |y| k )ρ(x, y, t) dx dy
Using integration by parts yields

dm k dt (t) = Dk(k -1)m k-2 (t) +k sgn(x)|x| k-1 F x [ρ] + sgn(y)|y| k-1 F y [ρ] ρ(x, y, t) dx dy = Dk(k -1)m k-2 (t) -k sgn(x)|x| k-1 sgn(x -x ′ )ρ(x ′ , y)ρ(x, y, t) dx dx ′ dy -k sgn(y)|y| k-1 sgn(y -y ′ )ρ(x ′ , y)ρ(x, y, t) dx dx ′ dy.
By exchanging the rôle of x and x ′ , we find

sgn(x)|x| k-1 sgn(x -x ′ )ρ(x ′ , y)ρ(x, y, t) dx dx ′ dy = 1 2 [sgn(x)|x| k-1 -sgn(x ′ )|x ′ | k-1 ]sgn(x -x ′ )ρ(x ′ , y, t)ρ(x, y, t) dx dx ′ dy ≥ 0 since x → sgn(x)|x| k-1 is non-decreasing.
A similar remark applies for the integral coming from F y . Therefore, the moments satisfy the following relation

dm k dt ≤ Dk(k -1)m k-2 . (7)
In particular, since the total mass is conserved

d dt ρ(x, y, t) dy dx = 0, we obtain m 2 (t) ≤ m 2 (0) + 2DM 0 t, with M 0 = ρ 0 (x, y) dy dx.

Entropies

Let h : [0, ∞) → R be a convex function and write

H[ρ] = h(ρ) dy dx. We have dH[ρ] dt = h ′ (ρ)∇ • D∇ρ -F [ρ]ρ dy dx = -D |∇ρ| 2 h ′′ (ρ) dy dx + ∇ρ • F [ρ] ρh ′′ (ρ) dy dx.
Let q be an anti-derivative of q ′ (ρ) = ρh ′′ (ρ).

We thus arrive at

dH[ρ] dt = -D |∇ρ| 2 h ′′ (ρ) dy dx - q(ρ)∇ • F [ρ] dy dx = -D |∇ρ| 2 h ′′ (ρ) dy dx + 4 ρq(ρ) dy dx (8)
by virtue of (4). In order to compensate the non-linearity in the last integral by the dissipated term, we can make use of the following Gagliardo-Nirenberg-Sobolev inequality (see e. g. [20, p. 125] or [5, Th. IX.9 with eq. ( 17) & eq. (85) p. 195]), which holds in R 2 for any p ≥ 1:

ξ p+1 dy dx ≤ C p ξ dy dx × |∇(ξ p/2 )| 2 dy dx. (9)
Let us detail how the estimates work in different cases:

• Entropy h(z) = z ln(z).

We get zq(z) = z 2 and we use [START_REF] Dalibard | Laser cooling of an optically thick gas: the simplest radiation pressure trap ?[END_REF] with p = 1. Remarking that |∇ρ| 2 ρ = 4|∇ √ ρ| 2 and taking into account the mass conservation, we are led to

d dt ρ ln(ρ) dy dx + 4(D -C 1 M 0 ) |∇ √ ρ| 2 dy dx ≤ 0. ( 10
)
It indicates a dissipation property when the diffusion coefficient is large enough

D ≥ C 1 M 0 .
Based on this, we can conjecture that solutions exist globally for large diffusion constants D.

• L p estimates: h(z) = z p . We get zq(z) = (p -1)z p+1 , and we use (9

) with p > 1. Remarking that h ′′ (ρ)|∇ρ| 2 = 4 p (p -1)|∇ρ p/2 | 2 , we are led to d dt ρ p dy dx + 4(p -1) D p -C p M 0 |∇ρ p/2 | 2 dy dx ≤ 0. (11)
Eq. [START_REF] Dobrušin | Vlasov equations[END_REF] shows that the L p norm of the solution is a non-increasing function of time when D is large enough, but how large depends on p with this approach. We are going to obtain finer estimates for large values of the constant p.

In order to eliminate the too restrictive condition on D, we use a different approach for the L p estimate. To this end, we use the Cauchy-Schwarz inequality and (9) and we obtain ρ p+1 dy dx ≤ ρ 2 dy dx From now on, we assume that

D > 2C 2 M 0 . (12)
Accordingly, the L 2 norm is dissipated and ρ 2 (t) dy dx ≤ ρ 2 0 dy dx holds. Therefore, we arrive at

d dt ρ p dy dx + 2D p -1 p |∇ρ p/2 | 2 dy dx ≤ K 1 p 2 ρ p dy dx (13)
with

K 1 = 2C 2 1 ρ0 2 L 2

D

. We use this relation to derive a L ∞ estimate, through an iterative argument on the exponent p which dates back to [START_REF] Alikakos | L p bounds of solutions of reaction-diffusion equations[END_REF]. Let us set

p k = 2 k , v k = ρ p k .
Let ω > 0. Eq. ( 13) tells us that

e -ωt d dt e ωt v k dy dx + 2D p k -1 p k |∇v k-1 | 2 dy dx ≤ (K 1 p 2 k + ω) v k dy dx.
We are going to estimate the right hand side by using the following Gagliardo-Nirenberg-Sobolev inequality (see e. g. [20, p. 125] or [5, eq. (85) p. 195])

ξ 2 dy dx ≤ C2 ξ dy dx |∇ξ| 2 dy dx 1/2 . ( 14
)
We combine this information with the Young inequality as follows

ξ 2 dy dx ≤ C2 ω 2 |∇ξ| 2 dy dx + C2 2ω ξ dy dx 2 .
We choose ω = ω k > 0 small enough to ensure

C2 (K 1 p 2 k + ω k ) ω k 2 ≤ D p k -1 p k . Since v k = v 2 k-1 , we are thus led to e -ω k t d dt e ω k t v k dy dx + D p k -1 p k |∇v k-1 | 2 dy dx ≤ C2 (K 1 p 2 k + ω k ) 2ω k v k-1 dy dx 2 .
Integrating from this relation, we obtain

v k (t) dy dx ≤ e -ω k t v k (0) dy dx + t 0 e -ω k s C2 (K 1 p 2 k + ω k ) 2ω k v k-1 (s) dy dx 2 ds ≤ e -ω k t v k (0) dy dx +e -ω k t e ω k t -1 ω k C2 (K 1 p 2 k + ω k ) 2ω k sup 0≤s≤t v k-1 (s) dy dx 2 .
In the right hand side, we make a convex combination appear, and we infer that

v k (t) dy dx ≤ max v k (0) dy dx, C2 (K 1 p 2 k + ω k ) 2ω 2 k sup 0≤s≤t v k-1 (s) dy dx 2 .

Let us set

L = max ρ 0 L 1 , ρ 0 L ∞ , δ k = C2 (K 1 p 2 k + ω k ) 2ω 2 k . Note that ω k behaves like 1 p 2 k
, and thus we can dominate δ k ≤ M p 6 k for some M > 0, so that, finally, we can find

A > 0 such that δ k ≤ A k . A direct recursion shows that v k (t) dy dx ≤ δ k δ p1 k-1 . . . δ p k-1 1 L p k which implies ρ(t) L p k ≤ L (A r k ) 1/p k , r k = k ℓ=0 (k -ℓ)p ℓ .
Since

r k p k = 1 2 k j=1 j 1 2 j-1 = 1 2 d ds 1 -s k+1 1 -s s=1/2 = 2 1 + (k + 2)e -(k+1) ln(2)
admits a finite limit as k → ∞, we deduce that the sequence ρ(t) L p k k∈N is bounded . The L ∞ bound follows by letting k go to ∞, and the bound depends on the initial L1 and L ∞ norms. The minimal D needed for this bound to be valid is unknown. We can recap our findings as follows.

Proposition 2.1 Let ρ be a sufficiently smooth solution of (1)-(3). Then, ρ satisfies the following properties:

i) mass is conserved ρ(t) dy dx = ρ 0 dy dx = M 0 , ii) if ρ 0 ∈ L p (R 2 ) and D > pC p M 0 , 1 then, ρ(t) L p ≤ ρ 0 L p , iii) if ρ 0 ∈ L 1 ∩ L ∞ (R 2
) and D > 2C 2 M 0 , then there exists a constant M > 0 such that 0 ≤ ρ(y, x, t) ≤ M holds for a.e. t ≥ 0, (x, y) ∈ R 2 . iv) if (x, y) → (x 2 + y 2 )ρ 0 (x, y) ∈ L 1 (R 2 ), then, for any t ≥ 0, (x, y) → (x 2 + y 2 )ρ(x, y, t) ∈ L 1 (R 2 ), and m 2 (t) ≤ m 2 (0) + 2DM 0 t.

Estimates for the regularized problem

To analyze the solutions of the regularized PDE ( 5)-( 6) and justify their convergence as ε tends to 0, we will need estimates uniform with respect to ε. The following proposition is the equivalent of Proposition 2.1 for the regularized solution.

Proposition 2.2 Let ρ (ε) ε>0 be the sequence of solutions of the regularized PDE (5)-( 6), associated to the initial data ρ (ε) 0 ε>0 . We assume that

ρ (ε) 0 ε>0 is bounded in L 1 (R 2 ) ∩ L ∞ (R 2 ).
Then, the following properties are satisfied:

i) mass is conserved ρ (ε) (t) dy dx = ρ (ε) 0 dy dx, ii) if D > pC p M 0 , then ρ (ε) is bounded in L ∞ (0, ∞; L p (R 2 )) and ρ (ε) (t) L p ≤ ρ (ε) 0 L p , iii) if D > 2C 2 M 0 , then ρ (ε) is bounded in L ∞ (0, ∞; L 2 (R 2 )) ∩ L ∞ ((0, ∞) × R 2 ), iv) if (x, y) → (x 2 + y 2 )ρ (ε) 0 (x, y) is bounded in L 1 (R 2 ), then (x 2 + y 2 )ρ (ε) (x, y, t) is bounded in L ∞ (0, T ; L 1 (R 2 )) for any 0 < T < ∞.
Proof. Item i) is clear. The proof of iv) repeats the same arguments as above, with a direct comparison to a pure diffusion. For ii) and iii), we will need the following consequence of the definition ( 6)

∂ x F (ε) x [ρ (ε) ](x, y, t) = -2 δ (ε) (x -x ′ )δ (ε) (y -y ′ )ρ (ε) (x -x ′ , y -y ′ , t) dx ′ dy ′ so that (4) becomes ∇ • F (ε) [ρ (ε) ] = -4T (ε) (ρ (ε) )
where, as said above, T (ε) is the convolution operator with the normalized 2d Gaussian kernel. Furthermore, the Hölder inequality yields

(ρ (ε) ) p T (ε) (ρ (ε) ) dy dx ≤ (ρ (ε) ) p+1 dy dx p/(p+1) |T (ε) (ρ (ε) )| p+1 dy dx 1/(p+1) ≤ (ρ (ε) ) p+1 dy dx p/(p+1) (ρ (ε) ) p+1 dy dx 1/(p+1) ≤ (ρ (ε) ) p+1 dy dx.
With this observation, we can go back to (8) adapted to the regularized problem and we derive the estimates as we did for the singular equation. We refer the reader to [START_REF] Blanchet | The two-dimensional Keller-Segel model: Optimal critical mass and qualitative properties of the solutions[END_REF] for similar reasonings.

Regularized problem

Let ε > 0. The initial data ρ

(ε) 0 is a given non-negative function in L 1 (R 2 ) ∩ L ∞ (R 2
). We introduce the operator

T : g → T (g) = ρ
where ρ is the solution of the linear parabolic PDE

∂ t ρ = D∆ρ -∇ • ( F (ε) [g]ρ), ρ t=0 = ρ (ε) 0 . ( 15
)
We will show that T fulfils the hypotheses of the Schauder theorem in a suitable functional framework. This will lead to the existence of a fixed point, which defines a solution of the non-linear problem. Then, we will investigate the uniqueness independently. Gathering together these arguments, we will prove the following statement.

Theorem 3.1 Let ρ (ε) 0 ∈ L 1 (R 2 ) ∩ L ∞ (R 2
) be a non-negative function. Then, the problem (5)-( 6)

with ρ t=0 = ρ (ε) 0 admits a unique solution ρ (ε) ∈ C([0, T ]; L 2 (R 2 )) ∩ L 2 (0, T ; H 1 (R 2 )).

Preparing for the Schauder theorem: a priori estimates

We observe that sgn

(ε) ∈ C ∞ (R) ∩ L ∞ (R), x → sgn (ε) (x) √ 1 + x 2 ∈ L 2 (R), δ (ε) = 1 2 d dx sgn (ε) ∈ C ∞ (R) ∩ L 1 (R) ∩ L ∞ (R).
Owing to these properties, we obtain estimates (that depend on ε) on the regularized force.

Lemma 3.2

The following estimates hold i)

|F (ε) x [g]| ≤ sgn (ε) ∞ δ (ε) ∞ g L 1 = 1 ε √ 2π g L 1 , ii) |F (ε) x [g]| ≤ √ π δ (ε) L 2 (1 + x ′2 )g 2 (x ′ , y ′ ) dx ′ dy ′ 1/2 , iii) |∂ x F (ε) x [g]| ≤ 2 δ (ε) 2 L 1 g L ∞ = 2 g L ∞ , iv) |∂ x F (ε) x [g]| ≤ 2 δ (ε) 2 L ∞ g L 1 = 2 πε 2 g L 1 .
Of course, the same estimates apply to F (ε) y as well.

Proof. It is worth bearing in mind that 0 ≤ δ (ε) (x) ≤ 1 ε √ 2π , |sgn (ε) (x)| ≤ 1.
Items i), iii) and iv) are direct consequences of estimates on convolution products. For ii) we use the Cauchy-Schwarz inequality twice to obtain

|F (ε) x [g](x, y)| ≤ δ (ε) (y -y ′ ) |sgn (ε) (x -x ′ )| 2 1 + x ′2 dx ′ 1/2 (1 + x ′2 )g 2 (x ′ , y ′ ) dx ′ 1/2 dy ′ ≤ |δ (ε) | 2 (y -y ′ ) dy ′ 1/2 |sgn (ε) (x -x ′ )| 2 1 + x ′2 dx ′ 1/2 (1 + x ′2 )g 2 (x ′ , y ′ ) dx ′ dy ′ 1/2 ≤ δ (ε) L 2 × √ π × (1 + x ′2 )g 2 (x ′ , y ′ ) dx ′ dy ′ 1/2 . For any g ∈ L ∞ (0, ∞; L 1 (R 2 )
), owing to the observations in Lemma 3.2, the linear problem (15) admits a unique solution, say in

C([0, ∞]; L 2 (R 2 )) ∩ L 2 (0, ∞; H 1 (R 2 )), see [5, Th. X.9].
We can now derive estimates on the solution of (15).

Lemma 3.3 Let ρ = T (g) be the solution of [START_REF] Hale | Asymptotic Behavior of Dissipative Systems[END_REF]. It satisfies i) For any fixed time t > 0 and any

p ∈ [1, ∞], ρ(t) ∈ L p (R 2 ). More precisely, we have ρ(x, y, t) p dy dx ≤ e 4(p-1)t g L ∞ (0,∞;L 1 (R 2 )) δ (ε) 2 L ∞ ρ (ε) 0 (x, y) p dy dx, and
ρ(t) L ∞ ≤ e 4t g L ∞ (0,∞;L 1 (R 2 )) δ (ε) 2 L ∞ ρ (ε) 0 L ∞ . ii)
For any fixed time t, (x 2 + y 2 )ρ(x, y, t) dy dx is finite. More precisely, we have

(x 2 + y 2 )ρ(x, y, t) dy dx ≤ e t (x 2 + y 2 )ρ (ε) 0 (x, y) dy dx + t ρ (ε) 0 L 1 4D + sgn (ε) 2 L ∞ δ (ε) 2 ∞ g 2 L ∞ (0,∞;L 1 (R 2 ))
.

iii) For any 0 ≤ t ≤ T < ∞, we have

∇ρ 2 L 2 ((0,t)×R 2 ) ≤ 1 2D e 4T g L ∞ (0,∞;L 1 (R 2 )) δ (ε) 2 L ∞ ρ (ε) 0 L 2 . Proof. i) We compute d dt ρ p dy dx + Dp(p -1) ρ p-2 |∇ρ| 2 dy dx = F (ε) [g]ρp(p -1) • ∇ρρ p-2 dy dx = -(p -1) ∇ • F (ε) [g]ρ p dy dx ≤ (p -1) ∇ • F (ε) [g] L ∞ ρ p dy dx ≤ 4 δ (ε) 2 L ∞ g L 1 × (p -1) ρ p dy dx.
The last line uses Lemma 3.2-iv). Grönwall's lemma then yields i). The L ∞ estimate follows by taking the limit p → ∞. Estimate iii) is obtained by specifying to the case p = 2 and considering the dissipation term.

ii. Let us use the shorthand notation z = (x, y). We get

d dt |z| 2 ρ dz = 4D ρ dz + 2 ρF (ε) [g] • z dz ≤ 4D ρ dz + 2 F (ε) [g] L ∞ |z| 2 ρ dz 1/2 ρ dz 1/2 ≤ 4D ρ dz + F (ε) [g] 2 L ∞ ρ dz + |z| 2 ρ dz ≤ (4D + sgn (ε) 2 ∞ δ (ε) 2 ∞ g 2 L 1 )
ρ dz + |z| 2 ρ dz by using Lemma 3.2-i). The Grönwall lemma allows us to conclude.

iii) We have

d dt ρ 2 dy dx + 2D |∇ρ| 2 dy dx ≤ 4 δ (ε) 2 L ∞ g L 1 × ρ 2 dy dx. ( 16
)
Inserting the estimate of item i) leads to

d dt ρ 2 dy dx + 2D |∇ρ| 2 dy dx ≤ Ce Ct ρ (ε) 0 (x, y) 2 dy dx, (17)
with

C = 4 δ (ε) 2 L ∞ g L ∞ (0,∞;L 1 (R 2 ))
. Integrating over time then yields

ρ 2 dy dx (t) - ρ 2 dy dx (t = 0) + 2D ∇ρ 2 L 2 ((0,t)×R 2 ) ≤ (e Ct -1) ρ (ε) 0 2 L 2 .( 18
)
Hence

∇ρ 2 L 2 ((0,t)×R 2 ) ≤ e Ct 2D ρ (ε) 0 2 L 2 (19)
≤ 1 2D e 4T g L ∞ (0,∞;L 1 (R 2 )) δ (ε) 2 L ∞ ρ (ε) 0 2 L 2 . ( 20
)

Preparing for the Schauder theorem: definition of the functional framework

Let 0 < T < ∞ be fixed once for all. We introduce the set C consisting of the functions g :

[0, T ] × R 2 → [0, ∞), such that i) g dz = ρ (ε) 0 dz ≤ M 0 , ii) |z| 2 g dz ≤ e C1T ( |z| 2 ρ (ε) 0 dz + C 2 T ), iii) g ∞ ≤ e C3T ρ (ε) 0 ∞ ,
By using the mass conservation property, the estimates in Lemma 3.3 allow us to choose the constants C 1 , C 2 and C 3 (which depend on ε) such that C is convex, and stable upon application of T .

Preparing for the Schauder theorem: T is continuous

We wish to establish the continuity of T :

C → C for the norm of L 2 ((0, T ) × R 2 ). For i ∈ {1, 2}, let ρ i = T (g i ), with g i ∈ C. By Lemma 3.3-i), we already know that ρ i belongs to L ∞ (0, T ; L 2 (R 2 )). We denote F (ε) i = F (ε) [g i ]. We get d dt (ρ 2 -ρ 1 ) 2 dz = -2D |∇(ρ 2 -ρ 1 )| 2 dz + 2 ∇(ρ 2 -ρ 1 ) • F (ε) 2 (ρ 2 -ρ 1 ) dz +2 ρ 1 ∇(ρ 2 -ρ 1 ) • ( F (ε) 2 -F (ε) 1 ) dz ≤ -2D |∇(ρ 2 -ρ 1 )| 2 dz + |∇ • F (ε) 2 | (ρ 2 -ρ 1 ) 2 dz +D |∇(ρ 2 -ρ 1 )| 2 dz + 1 D ρ 2 1 ( F (ε) 2 -F (ε) 1 ) 2 dz ≤ -D |∇(ρ 2 -ρ 1 )| 2 dz + 4 δ (ε) 2 L ∞ g 2 L 1 (ρ 2 -ρ 1 ) 2 dz + 1 D F (ε) 2 -F (ε) 1 2 L ∞ ρ 2 1 dz. ( 21
)
We aim at controlling

F (ε) 2 -F (ε) 1 2
L ∞ by the difference g 2 -g 1 in L 2 norm. This cannot be done directly, and we should use further moment estimates. To be more specific, we will use a splitting that makes g 2 -g 1 L 2 appear plus an arbitrarily small contribution. To this end, we use Lemma 3.2-ii). For any 1 < s < 2 and any R > 0, we write (F

(ε) 2,x -F (ε) 1,x ) 2 (x, y, t) ≤ δ (ε) 2 L 2 |sgn (ε) (x -x ′ )| 2 1 + |x ′ | s dx ′ (1 + |x ′ | s )(g 2 -g 1 ) 2 (x ′ , y ′ , t) dx ′ dy ′ ≤ C (ε) (1 + |z| s )(g 2 -g 1 ) 2 (z, t) dz ≤ C (ε) (g 2 -g 1 ) 2 (z, t) dz + |z|≤R |z| s (g 2 -g 1 ) 2 (z, t) dz + |z|>R |z| s (g 2 -g 1 ) 2 (z, t) dz ≤ C (ε) (1 + R s ) (g 2 -g 1 ) 2 (z, t) dz + C (ε) |z|>R |z| 2 |z| 2-s (g 2 -g 1 ) 2 (z, t) dz ≤ C (ε) (1 + R s ) (g 2 -g 1 )(t) 2 L 2 + C (ε) g 2 L ∞ + g 1 L ∞ R 2-s |z| 2 |g 2 -g 1 |(z, t) dz ≤ C (ε) (1 + R s ) (g 2 -g 1 )(t) 2 L 2 + C (ε) g 2 L ∞ + g 1 L ∞ R 2-s |z| 2 g 2 dz + |z| 2 g 1 dz .
The same inequalities obviously hold for

F (ε) 2,y -F (ε) 1,y . Coming back to (21) yields d dt (ρ 2 -ρ 1 ) 2 dz + D |∇(ρ 2 -ρ 1 )| 2 dz ≤ 4 δ (ε) 2 L ∞ g 2 L 1 (ρ 2 -ρ 1 ) 2 dz + 2C (ε) ρ 1 (t) 2 L 2 D (1 + R s ) (g 2 -g 1 )(t) 2 L 2 + g 2 L ∞ + g 1 L ∞ R 2-s |z| 2 g 2 + |z| 2 g 1 .
Bearing in mind that C 3 = 4 δ (ε) 2 ∞ M 0 , we are ready to use the Grönwall lemma which leads to

(ρ 1 -ρ 2 ) 2 (z, t) dz ≤ e C3T (ρ 2 -ρ 1 ) 2 (z, 0) dz + 2C (ε) ρ 1 2 L ∞ (0,T ;L 2 (R 2 )) (1 + R s ) D t 0 (g 2 -g 1 ) 2 (z, τ ) dz dτ + 2C (ε) ρ 1 2 L ∞ (0,T ;L 2 (R 2 )) ( g 2 L ∞ + g 1 L ∞ ) DR 2-s t 0 |z| 2 (g 2 + g 1 )(z, τ ) dz dτ .
When ρ 1 and ρ 2 have the same initial condition the first term of the right hand side vanishes. Take g ∈ C and consider a sequence

g n n∈N ∈ C, such that g n → g in L 2 ([0, T ] × R 2
). We apply [START_REF] Poupaud | Parabolic limit and stability of the Vlasov-Poisson-Fokker-Planck system[END_REF] with ρ n = T (g n ) and ρ = T (g); it reads

(ρ n -ρ) 2 (z, t) dz ≤ e C3T 2C (ε) ρ 2 L ∞ (0,T ;L 2 (R 2 )) (1 + R s ) D T 0 (g n -g) 2 (z, τ ) dz dτ +e C3T 2C (ε) ρ 2 L ∞ (0,T ;L 2 (R 2 )) ( g n L ∞ + g L ∞ ) DR 2-s T 0 |z| 2 (g n + g)(z, τ ) dz dτ
Pick η > 0. Using the bounds that define the set C, it is possible to select R(η) > 0 such that the last term can be made smaller than η/2, uniformly with respect to n. Then, with this R at hand, there exists N (η) ∈ N such that for all n ≥ N (η) the first term in the right hand side is smaller than η/2 too. Hence,

(ρ n -ρ) 2 (z, t) dz ≤ η holds for any n ≥ N (η) and 0 ≤ t ≤ T < ∞. It shows that ρ n → ρ in L 2 ((0, T ) × R 2 ). Thus T : C → C is continuous for the strong topology of L 2 ((0, T ) × R 2 ).

Preparing for the Schauder theorem: T is compact

Let g n n∈N be a sequence in C, and set

ρ n = T (g n ). Then, from Lemma 3.3, ρ n is bounded in L ∞ 0, T ; L 2 (R 2 )
, and, furthermore, ∇ρ n is bounded in L 2 (0, T ) × R 2 . We also have

∂ t ρ n = D∇ • ∇ρ n -∇ • F (ε) [g n ]ρ n
where, by Lemma 3.2-i),

F (ε) [g n ] is bounded in L ∞ uniformly with respect to n. Therefore, ∂ t ρ n is bounded in L 2 0, T ; H -1 (R 2 ) . Since the embedding H 1 (B(0, R)) ⊂ L 2 (B(0, R)) is compact for any 0 < R < ∞,
|ρ n | 2 dz dt ≤ ρ n ∞ R 2 T 0 |z| 2 ρ n dz dt ≤ C(ε, T ) R 2
can be made arbitrarily small by choosing R large enough, uniformly with respect to n ∈ N. The sequence ρ n n∈N thus fulfils the criterion of the Fréchet-Weil-Kolmogorov theorem, see e. g [14, Th. [START_REF] Chalony | Long-range one-dimensional gravitational-like interaction in a neutral atomic cold gas[END_REF].56] and it is thus relatively compact in L 2 ((0, T ) × R 2 ).

Schauder theorem: existence

Gathering the results of the previous subsections, we can use Schauder's theorem:

C is a closed convex subset of L 2 ((0, T ) × R 2 ), T is a continuous mapping such that T (C) ⊂ C and T (C) is relatively compact in L 2 ((0, T )× R 2 )
. Then T has a fixed point, which is a solution of the nonlinear regularized PDE ( 5)-( 6) with initial condition ρ

(ε) 0 , on any arbitrary time interval [0, T ]. The obtained solution lies in C([0, T ]; L 2 (R 2 )) ∩ L 2 (0, T ; H 1 (R 2 )).

Uniqueness

The argument to justify uniqueness relies on the following claim, for which we refer the reader to [15, Lemma 7.1.1] or [START_REF] Dixon | Weakly singular discrete Gronwall inequalities[END_REF]Th. 3.1].

Lemma 3.4 (Singular Grönwall Lemma) Let A, B ≥ 0, 0 ≤ α < 1. Let u(t) a locally bounded function such that u(t) ≤ A + B t 0 u(s) (t -s) α ds then we have u(t) ≤ AE 1-α BΓ(1 -α)t 1-α
with s → Γ(s) the usual Γ-function and E 1-α stands for the Mittag-Leffler function with parameter

β = 1 -α E β = ∞ n=0 s n Γ(nβ + 1)
.

Let ρ 1 and ρ 2 be two solutions of the regularized nonlinear PDE. Let

H t (z) = 1 4πDt e -|z| 2 /(4Dt) (22)
stand for the two-dimensional heat kernel with coefficient D. We write

(ρ 1 -ρ 2 )(t) = H t ⋆ (ρ 1 -ρ 2 )(0) - t 0 H t-s ⋆ ∇ • ( F (ε) [ρ 2 ]ρ 2 -F (ε) [ρ 1 ]ρ 1 )(s) ds = H t ⋆ (ρ 1 -ρ 2 )(0) + t 0 ∇H t-s ⋆ ( F (ε) [ρ 2 ]ρ 2 -F (ε) [ρ 1 ]ρ 1 )(s) ds.
Initially we have ρ 2 (0) = ρ 1 (0) and (with

C 0 = 1 π |z|e -|z| 2 dz) we arrive at |ρ 1 -ρ 2 |(z, t) dz ≤ t 0 |∇H t-s (z -z ′ )| | F (ε) [ρ 2 ]ρ 2 -F (ε) [ρ 1 ]ρ 1 |(z ′ , s) dz ′ dz ds ≤ t 0 C 0 √ t -s | F (ε) [ρ 2 ]| |(ρ 2 -ρ 1 ))|(z ′ , s) dz ′ ds + t 0 C 0 √ t -s | F (ε) [ρ 2 ] -F (ε) [ρ 1 ]| ρ 1 (z ′ , s)| dz ′ ds. ( 23
)
Lemma 3.2-i) together with the mass conservation tell us that

F (ε) [ρ 2 ] L ∞ ≤ 2M 0 δ (ε) L ∞ sgn (ε) L ∞ .
We also have

ρ 1 (z, t)| F (ε) [ρ 2 ] -F (ε) [ρ 1 ]|(z, t) dz ≤ 2 δ (ε) L ∞ sgn (ε) L ∞ ρ 1 (z, t) dz |ρ 1 -ρ 2 |(z ′ , t) dz ′ ≤ 2M 0 δ (ε) L ∞ sgn (ε) L ∞ |ρ 1 -ρ 2 |(z ′ , t) dz ′ .
Introducing this into (23) yields, for a certain constant B > 0:

ρ 1 -ρ 2 L 1 (t) ≤ B t 0 1 √ t -s ρ 1 -ρ 2 L 1 (s)ds (24)
The singular Grönwall lemma allows us to conclude that

ρ 1 = ρ 2 .
4 Convergence of ρ (ε)

We can now state our main result about the existence of solutions for ( 1)-( 3), which is expressed as a stability result.

Theorem 4.1 Let ρ (ε)
0 be a sequence of non negative functions bounded in L 1 (R 2 ) and in L ∞ (R 2 ). We suppose that ρ (ε) 0 dy dx ≤ M 0 and D > 2C 2 M 0 (see [START_REF] Fournier | On the rate of convergence in Wasserstein distance of the empirical measure[END_REF]). Then, up to a subsequence, the associated sequence ρ (ε) ε>0 converges strongly in L p ((0, T ) × R 2 ) for any 1 ≤ p < ∞, and in C([0, T ]; L p (R 2 ) -weak), to ρ, which is a solution of (1)-( 2) with initial data ρ 0 , the weak limit of ρ (ε) 0 .

Compactness approach

We remind the reader that we are assuming D > 2C 2 M 0 . Accordingly, from Proposition 2.2, we already know that ρ (ε) ε>0 is bounded in L ∞ ((0, T ); L p (R 2 )), for any 1 ≤ p ≤ ∞, and ∇ρ (ε) is bounded in L 2 ((0, T ) × R 2 ). Moreover, the equation

∂ t ρ (ε) = ∇ • D∇ρ (ε) -F (ε) [ρ (ε) ]ρ (ε)
tells us that ∂ t ρ (ε) is the space derivative of the sum of a term bounded in L 2 ((0, T ) × R 2 )) and the divergence of a term bounded in L ∞ (0, T ; L 1 (B(0, R))) for any 0 < R < ∞. Indeed, we readily check that

B(0,R) |F (ε) x [ρ (ε) ](x, y, t)| dy dx ≤ √ x 2 +y 2 ≤R δ (ε) (y -y ′ )ρ (ε) (x ′ , y ′ , t) dx ′ dy ′ dy dx ≤ +R -R dx δ (ε) (y -y ′ ) dy ρ (ε) (x ′ , y ′ , t) dx ′ dy ′ (25) ≤ 2RM 0 . (26)
In fact it turns out that

F (ε) x [ρ (ε) ], like F x [ρ], is bounded in L ∞ ((0, T ) × R; L 1 (R)). Hence, ∂ t ρ (ε)
is bounded in, say, L 2 (0, T ; H -1-δ (B(0, R))) for any 0 < R < ∞ and δ > 0. We can apply the Aubin-Simon lemma [START_REF] Simon | Compact sets in L p (0, T ; B). Ann. Mat. Pura e Applicada[END_REF] and we conclude that ρ (ε) ε>0 is relatively compact in L 2 ((0, T ) × B(0, R)) for any 0 < T, R < ∞. By using the moments estimate, and reasoning as we did in Section 3.4, we show that ρ (ε) is actually relatively compact in L 2 ((0, T ) × R 2 ).

Therefore, possibly at the price of extracting a subsequence (still labelled by ε, though) we can assume that ρ (ε) → ρ strongly in L 2 ((0, T ) × R 2 ).

The convergence can be strengthened in two directions. First of all, if 1 < p = θ2 + (1 -θ) < 2, the Hölder inequality leads to

ρ (ε) -ρ L p ((0,T )×R 2 ) ≤ (2M 0 ) 1-θ ρ (ε) -ρ θ L 2 ((0,T )×R 2 ) and if 2 < p < ∞, we have ρ (ε) -ρ L p ((0,T )×R 2 ) ≤ ( ρ (ε) L ∞ + ρ L ∞ ) (p-2)/p ρ (ε) -ρ 2/p L 2 ((0,T )×R 2 )
. We can also treat the case p = 1 since the L 2 estimate and the moment estimate imply that ρ (ε) ε>0 is weakly compact in L 1 ((0, T ) × R 2 ) and we can assume that it converges a.e., see [START_REF] Goudon | Intégration. Intégrale de Lebesgue et introduction à l'analyse fonctionnelle[END_REF]Th. 7.60]. Finally we get

ρ (ε) → ρ strongly in L p ((0, T ) × R 2 ) for any 1 ≤ p < ∞. ( 27
)
Second of all, the bound on ∂ t ρ (ε) can be used to justify, by using the Arzela-Ascoli theorem and a diagonal extraction, that lim ε→0 ρ (ε) (x, y, t)φ(x, y) dy dx = ρ(x, y, t)φ(x, y) dy dx holds for any φ ∈ C(R 2 ), or in L p ′ (R 2 ), uniformly on [0, T ]. In particular, the initial data passes to the limit and (3) makes sense (with ρ 0 the weak limit in L p (R 2 ) of the extracted sequence ρ (ε) 0 ε>0 ). We are left with the task of passing to the limit in the non-linear term F (ε) [ρ (ε) ]ρ (ε) . To this end, we split as follows

F (ε) [ρ (ε) ] -F [ρ] = F (ε) [ρ (ε) -ρ] + ( F (ε) [ρ] -F [ρ]).
The first term tends to 0 as a consequence of ( 27) combined with the following claim.

Lemma 4.2 The operator

F (ε) x (resp. F (ε) y ) is, uniformly with respect to ε, continuous from L 1 (R 2 ) to L ∞ (R x ; L 1 (R y )) (resp. L ∞ (R y ; L 1 (R x ))).
Proof. For any φ ∈ L 1 (R 2 ), we have

|F (ε) x [φ](x, y)| dy ≤ δ (ε) (y -y ′ )|φ(x ′ , y ′ )| dx ′ dy ′ dy ≤ δ (ε) (y -y ′ ) dy |φ(x ′ , y ′ )| dx ′ dy ′ = φ L 1 .
It remains to investigate, for φ ∈ L 1 (R 2 ), the behavior of

F (ε) x [φ] -F x [φ] (x, y) = δ (ε) (y -y ′ )sgn (ε) (x -x ′ )φ(x ′ , y ′ ) dy ′ dx ′ - δ (ε) (y -y ′ ) dy ′ sgn(x -x ′ )φ(x ′ , y) dx ′ ≤ δ (ε) (y -y ′ ) sgn (ε) (x -x ′ )φ(x ′ , y ′ ) -sgn(x -x ′ )φ(x ′ , y) dy ′ dx ′ .
We integrate with respect to y and, bearing in mind that δ (ε) (y) = 1 ε δ(y/ε) with δ the normalized Gaussian, we use the change of variable y -y ′ = εξ; it yields

F (ε) x [φ] -F x [φ] (x, y) dy ≤ δ(ξ) sgn (ε) (x -x ′ )φ(x ′ , y -εξ) -sgn(x -x ′ )φ(x ′ , y) dξ dx ′ dy ≤ δ(ξ) φ(x ′ , y -εξ) -φ(x ′ , y) dξ dx ′ dy + δ(ξ)φ(x ′ , y) sgn (ε) (x -x ′ ) -sgn(x -x ′ ) dξ dx ′ dy.
On the right hand side, the first integral recasts as δ(ξ) φ(x ′ , y -εξ) -φ(x ′ , y) dx ′ dy, dξ which tends to 0 as ε → 0 by combining the Lebesgue dominated convergence theorem with the continuity of translation in L 1 , [START_REF] Goudon | Intégration. Intégrale de Lebesgue et introduction à l'analyse fonctionnelle[END_REF]Cor. 4.14]. The second integral reads

δ(ξ) dξ × φ(x ′ , y) sgn (ε) (x -x ′ ) -sgn(x -x ′ ) dx ′ dy The function (x, x ′ ) → |sgn (ε) (x -x ′ ) -sgn(x -x ′ )
| tends to 0 pointwise and it is dominated by 2. Since ρ ∈ L 1 (R 2 ), a direct application of the Lebesgue dominated convergence theorem tells us that this quantity tends to 0 as ε → 0, for any given x ∈ R. Similar reasoning obviously apply to the second component of F . Finally, for any test function

ϕ ∈ C ∞ c (R 2 ), we obtain lim ε→0 ϕ F (ε) [ρ (ε) ]ρ (ε) -F [ρ]ρ dy dx = 0.
Therefore ρ satisfies, in a weak sense, the limit equation ( 1)-( 2). This completes the proof of Theorem 4.1.

Symmetric solutions

Throughout this Section, we work with data that satisfy the following symmetry condition ρ 0 (-x, y) = ρ 0 (x, y) = ρ 0 (x, -y).

(28)

It will be used to derive further estimates and a stronger convergence result of the regularized solution ρ (ε) towards the solutions of (1)- [START_REF] Blanchet | The two-dimensional Keller-Segel model: Optimal critical mass and qualitative properties of the solutions[END_REF]. Using the uniqueness property of the solution of the regularized equation ( 5)-( 6), we deduce that the symmetry property is preserved by the solutions of (1). Accordingly, we get

F x [ρ](0, y, t) = -sgn(x ′ )ρ(x ′ , y, t) dx ′ = 0, F y [ρ](x, 0, t) = 0.
However, we know that

∂ x F x [ρ] < 0 and ∂ y F y [ρ] < 0. Thus, x → F x [ρ](x, y, t)
is non increasing and it vanishes for x = 0, so that it has the sign of (-x). We deduce that

(x, y) • F [ρ](x, y, t) = xF x [ρ](x, y, t) + yF y [ρ](x, y, t) ≤ 0.
A similar property hold with the solutions ρ (ε) of the regularized problem and the force operator ε) ]. This will be used to obtain a strengthened control on the behavior of the solutions for large x, y's: exponential moments and weighted estimates on the gradients. These estimates will be combined with the interpretation of (1) as a perturbation of the heat equation. Namely, still with H t the heat kernel [START_REF] Poupaud | Parabolic limit and stability of the Vlasov-Poisson-Fokker-Planck system[END_REF], we shall make use of the Duhamel formula

F (ε) [ρ (
ρ(x, y, t) = H t ⋆ ρ 0 (x, y) - t 0 H t-s ⋆ ∇ • F [ρ]ρ(s, •) (x, y) ds, (29)
and the analogous formula with ρ (ε) .

Strengthened estimates for symmetric solutions

At first, the symmetry property allows us to control exponential moments.

Lemma 4.3 (Exponential moments) Assume that ρ 0 satisfies (28) and e λ √ 1+x 2 +y 2 ρ 0 (x, y) dy dx = E 0 (λ) < +∞ for some λ > 0. Then, the solutions of (1)-( 3) satisfy e λ √ 1+x 2 +y 2 ρ(x, y, t) dy dx ≤ E 0 (λ)e D(λ 2 +2λ)t .

The same estimate holds replacing ρ by ρ (ε) .

Proof As consequence of the symmetry assumption, the last term contributes negatively. We end the proof by integrating with respect to time.

Using L q and moments estimates, we can readily obtain a weighted L 2 bound; for instance, we have where the constant C depends on D, E 0 (2λ), ρ 0 L1 and ρ 0 L ∞ . Again, the same (uniform) estimate is fulfilled by ρ (ε) .

e λ √ 1+x
We need now to specify the class of initial data to which the analysis applies. Addtionally to the symmetry assumption, we suppose that ρ (ε) 0 , which is a regularization of ρ 0 in (3), is such that there exists p 0 , p 2 ≥ 0 such that for any λ ≥ 0, we have

E 0 (λ) = sup ε>0 e λ √ 1+x 2 +y 2 ρ (ε) 0 (x, y) dy dx ≤ e p0+p2λ 2 . (31)
Such an assumption clearly holds for uniformly compactly supported data, as well as for Gaussianlike data. Finally, for our purpose, we will need another estimate for the weighted L 2 norm of the gradient, which applies for the data verifying (31).

0 dz + (1 + |z| 2 )|ρ (ε) 0 | 2 dz < ∞.
If D satisfies (12), there exists constants B 0 , B 1 , B 2 > 0, which do not depend on ε nor on t, such that

(1 + |z| 2 )(ρ (ε) ) 2 dz ≤ B 0 + B 1 t + B 2 t 2 , t 0 (1 + |z| 2 )|∇ρ (ε) (z, s)| 2 ds ≤ B 0 + B 1 t + B 2 t 2 .
Proof. Let us compute (still with the shorthand notation z = (x, y)), by using several integrations by parts,

1 2 d dt (1 + |z| 2 )|ρ (ε) | 2 (z, t) dz = -D (1 + |z| 2 )|∇ρ (ε) | 2 dz + 2D |ρ (ε) | 2 dz + |ρ (ε) | 2 z • F (ε) [ρ (ε) ] dz - 1 2 (1 + |z| 2 )|ρ (ε) | 2 ∇ • ( F (ε) ρ (ε) ) dz.
The symmetry assumption implies z • F (ε) [ρ (ε) ] ≤ 0, which allows us to get rid of the third term in the right side. We remind the reader that ∇ • (

F (ε) ρ (ε) ) = -4T (ε) (ρ (ε)
) is proportional to the convolution with an approximation of the 2D Dirac measure. Hence, we get

1 2 d dt (1 + |z| 2 )|ρ (ε) | 2 (z, t) dz + D (1 + |z| 2 )|∇ρ (ε) | 2 dz ≤ 2D |ρ (ε) | 2 dz + 2 (1 + |z| 2 )|ρ (ε) | 2 T (ε) (ρ (ε) ) dz ≤ 2D |ρ (ε) | 2 dz + 2 (1 + |z| 2 ) 2 ρ (ε) dz 1/2 T (ε) (ρ (ε) ) 2 ρ (ε) 3 dz 1/2 ≤ 2D |ρ (ε) | 2 dz + 2 2 (1 + |z| 4 )ρ (ε) dz 1/2 |ρ (ε) | 5 dz 1/2
.

For the last term, we have used Hölder's inequality as in the proof of Proposition 2.2. We already know that the L 2 and L 5 norms of ρ (ε) are uniformly bounded, by virtue of Proposition 2.2. It remains to discuss the forth order moment. To this end we go back to (7): t → m 2 (t) has a linear growth, hence t → m 4 (t) have a quadratic growth with respect to the time variable. We conclude that both

(1 + |z| 2 )|ρ (ε) | 2 (z, t) dz and t 0 (1 + |z| 2 )|∇ρ (ε) (z, s)| 2 dz ds
has at most a quadratic growth, with coefficients independent of ε.

Cauchy property for ρ (ε)

This Section is concerned with the following statement, which strengthens Theorem 4.1 for symmetric solutions.

Theorem 4.5 Let ρ (ε) 0 ε>0 be a sequence of non negative functions bounded in L 1 ∩ L ∞ (R 2 ), which satisfies (28) and (31) and which converges in L 1 (R 2 ) to some ρ 0 . Then the associated sequence ρ (ε) ε>0 of solutions of (5)-( 6) is a Cauchy sequence in C([0, T ]; L 1 (R 2 )).

Corollary 4.6 Let ρ 0 ∈ L 1 ∩L ∞ (R 2 ) verify (28) and (31). Then the sequence ρ (ε) ε>0 of solutions of (5)-( 6) with the same initia data converges C([0, T ]; L 1 (R 2 )) to ρ, the unique symmetric solution of (1)-(3).

We make use of (29), which leads to

|ρ (ε) -ρ (ε ′ ) |(z, t) dz ≤ H t ⋆ |ρ (ε) 0 -ρ (ε ′ ) 0 |(z) dz + t 0 ∇H t-s ⋆ F (ε) [ρ (ε) ]ρ (ε) -F (ε ′ ) [ρ (ε ′ ) ]ρ (ε ′ ) (s, •) (z) dz ds.
(32) We dominate the right hand side by the sum of the following four terms

A ε,ε ′ (t) = R 2 H t ⋆ |ρ (ε) 0 -ρ (ε ′ ) 0 |(z) dz, B ε,ε ′ (t) = t 0 |∇H t-s (z -z ′ )| | F (ε) [ρ (ε) ](z ′ , s)| |ρ (ε) -ρ (ε ′ ) |(z ′ , s) dz ′ dz ds, C ε,ε ′ (t) = t 0 |∇H t-s (z -z ′ )| ρ (ε ′ ) (z ′ , s) | F (ε) [ρ (ε) -ρ (ε ′ ) ](z ′ , s)| dz ′ dz ds, D ε,ε ′ (t) = t 0 ρ (ε ′ ) (z ′ , s) ∇H t-s (z -z ′ ) • F (ε) [ρ (ε ′ ) ] -F (ε ′ ) [ρ (ε ′ ) ](z ′ , s) dz ′ dz ds. Since ρ (ε) 0 → ρ 0 in L 1 (R 2 ), it is clear that lim ε,ε ′ →0 sup t≥0 A ε,ε ′ (t) = 0 (33)
uniformly on any time interval [0, T ]. Next, we are going to justify the following claim.

Lemma 4.7 Let α = 1 √ 2 . Set ϕ(λ) = 2Dλ(1 + α) 1 + λ(1 + α) .
Then there exists constant β 1 , β 2 > 0 such that, for any R > 0 we have

B ε,ε ′ (t) ≤ β 1 R ρ (ε) L ∞ t 0 1 √ t -s (ρ (ε) -ρ (ε ′ ) )(s, •) L 1 ds + β 2 √ t λ e -αλR e ϕ(λ)t . ( 34
)
The constant β 1 does not depend on the data, while β 2 depends on E 0 (2λ(1 + α)).

Proof. In Section 3.6, we already used the basic estimate

|∇H t-s (z -z ′ )| dz ′ ≤ C 0 √ t -s (35)
for a certain constant C 0 . We have

|F (ε) [ρ (ε) ](z ′ )||ρ (ε) -ρ (ε ′ ) |(z ′ ) dz ′ ≤ ρ (ε) (x 1 , y 1 )δ (ε) (y ′ -y 1 )|ρ (ε) -ρ (ε ′ ) |(x ′ , y ′ ) dx ′ dy ′ dx 1 dy 1 ≤ x 2 1 +y 2 1 ≤R 2 ... dx ′ dy ′ dx 1 dy 1 + x 2 1 +y 2 1 ≥R 2 ... dx ′ dy ′ dx 1 dy 1 .
We dominate the first integral as follows

x 2 1 +y 2 1 ≤R 2 ... dx ′ dy ′ dx 1 dy 1 ≤ ρ (ε) L ∞ x1≤R δ (ε) (y ′ -y 1 ) dy 1 |ρ (ε) -ρ (ε ′ ) |(x ′ , y ′ ) dy ′ dx ′ ≤ 2R ρ (ε) L ∞ ρ (ε) -ρ (ε ′ ) L 1 .
Next, we have

x 2 1 +y 2 1 ≥R 2 ... dx ′ dy ′ dx 1 dy 1 ≤ x 2 1 +y 2 1 ≥R 2 ρ (ε) (x 1 , y 1 )ρ (ε) (x ′ , y ′ )δ (ε) (y ′ -y 1 ) dx ′ dy ′ dx 1 dy 1 + x 2 1 +y 2 1 ≥R 2 ρ (ε) (x 1 , y 1 )ρ (ε ′ ) (x ′ , y ′ )δ (ε) (y ′ -y 1 ) dx ′ dy ′ dx 1 dy 1
where the two terms can be treated with the same approach. We make the exponential moment appear and we use the Cauchy-Schwarz inequality to obtain, for instance,

x 2 1 +y 2 1 >R 2 ρ (ε) (x 1 , y 1 )δ (ε) (y -y 1 )ρ (ε) (x, y) dx 1 dy 1 dx dy ≤ 1 2 x 2 1 +y 2 1 >R 2 e λ √ 1+x 2 1 +y 2 1 e -λ √ 1+x 2 +y 2 |ρ (ε) | 2 (x 1 , y 1 )δ (ε) (y -y 1 ) dx 1 dy 1 dx dy + 1 2 x 2 1 +y 2 1 >R 2 e -λ √ 1+x 2 1 +y 2 1 e λ √ 1+x 2 +y 2 |ρ (ε) | 2 (x, y)δ (ε) (y -y 1 ) dx 1 dy 1 dx dy.
(36) The elementary inequality

α(|x| + |y|) ≤ 1 + x 2 + y 2 allows us to estimate e -λ √ 1+x 2 +y 2 dx ≤ e -αλ|y| 2 αλ .
Hence the first integral in the right hand side of (36) is dominated by

x 2 1 +y 2 1 >R 2 δ (ε) (y -y 1 ) e -λ √ 1+x 2 +y 2 dx dy e λ √ 1+x 2 1 +y 2 1 |ρ (ε) | 2 (x 1 , y 1 ) dx 1 dy 1 ≤ 2 αλ x 2 1 +y 2 1 >R 2 δ (ε) (y -y 1 )e -αλ|y| dy e λ √ 1+x 2 1 +y 2 1 |ρ (ε) | 2 (x 1 , y 1 ) dx 1 dy 1 ≤ 2 αλ x 2 1 +y 2 1 >R 2 e λ √ 1+x 2 1 +y 2 1 |ρ (ε) | 2 (x 1 , y 1 ) dx 1 dy 1 ≤ 2 αλ x 2 1 +y 2 1 >R 2 e -λα √ 1+x 2 1 +y 2 1 e λ(1+α) √ 1+x 2 1 +y 2 1 |ρ (ε) | 2 (x 1 , y 1 ) dx 1 dy 1 ≤ 2 αλ e -λαR e λ(1+α) √ 1+x 2 1 +y 2 1 |ρ (ε) | 2 (x 1 , y 1 ) dx 1 dy 1 ≤ C 2 αλ e -λαR e 2Dλ(1+α)(1+(1+α)λ)t
where we have used (30) and the constant C > 0 depends on E 0 (2λ(1 + α)). Next, we observe that

1 x 2 1 +y 2 1 >R 2 e -λ √ 1+x 2 1 +y 2 1 dx 1 ≤ e -αλ|y1| 1 x 2 1 +y 2 1 >R 2 e -αλ|x1| dx 1 ≤ 2 αλ e -αλ(|y1|+g(y1))
where

g(y) = 1 |y|≤R R 2 -y 2 .
As a matter of fact, for any y ∈ R, we have |y| + g(y) ≥ R, so that

1 x 2 1 +y 2 1 >R 2 e -λ √ 1+x 2 1 +y 2 1 dx 1 ≤ 2 αλ e -αλR .
The second integral of the right hand side in (36), is thus dominated by

1 x 2 1 +y 2 1 >R 2 e -λ √ 1+x 2 1 +y 2 1 dx 1 δ (ε) (y -y 1 ) dy 1 e λ √ 1+x 2 +y 2 |ρ (ε) | 2 (x, y) dx dy ≤ 2 αλ e -αλR δ (ε) (y -y 1 ) dy 1 e λ √ 1+x 2 +y 2 |ρ (ε) | 2 (x, y) dx dy ≤ C ′ 2 αλ e -αλR e 2Dλ(1+λ)t
where we have used (30) again and C ′ here depends on E 0 (2λ). We finally conclude (note that λ(1 + α) > λ) that

|F (ε) [ρ (ε) ](z ′ , s)||ρ (ε) -ρ (ε ′ ) |(z ′ , s) dz ′ ≤ 2R ρ (ε) L ∞ (ρ (ε) -ρ (ε ′ ) )(s) L 1 + C αλ e -αλR e ϕ(λ)t
with C depending on E 0 (2λ(1 + α)) ≥ E 0 (2λ). We combine this inequality to (35) to obtain the final estimate on B ε,ε ′ .

Lemma 4.8 There exists constant γ 1 , γ 2 > 0 such that, for any R > 0 we have

C ε,ε ′ (t) ≤ γ 1 R ρ (ε) L ∞ t 0 1 √ t -s (ρ (ε) -ρ (ε ′ ) )(s, •) L 1 ds + γ 2 √ t λ e -αλR e ϕ(λ)t . ( 37
)
The constant γ 1 does not depend on the data, while γ 2 depends on E 0 (2λ(1 + α)).

Proof. The same reasoning applies for C ε,ε ′ . Indeed, we can first integrate ∇H t-s (z -z ′ ) over z, which leads to the analog of (35). Estimating F (ε) , we are left with

C ε,ε ′ (t) ≤ t 0 C 0 √ t -s ρ (ε ′ ) (x ′ , y ′ , s)|ρ (ε) -ρ (ε ′ ) |(x 1 , y 1 , s))δ (ε) (y ′ -y 1 ) dx 1 dy 1 dx ′ dy ′ ds
which is exactly the same expression that appeared in the analysis of B ε,ε ′ .

We turn to the analysis of D ε,ε ′ . Lemma 4.9 Let 0 < T < ∞. Then D ε,ε ′ (t) converges to 0, uniformly over [0, T ] as ε, ε ′ tend to 0

Proof. We evaluate D ε,ε ′ through the following splitting

D ε,ε ′ (t) ≤ D x,1 (t) + D x,2 (t) + D y,1 (t) + D y,2 (t) with D x,1 (t) = t 0 ∂ x H t-s (x -x ′ , y -y ′ )ρ (ε ′ ) (x ′ , y ′ , s)ρ (ε ′ ) (x", y", s) × sgn (ε) (x ′ -x") δ (ε) (y ′ -y") -δ (ε ′ ) (y ′ -y") dx ′ dy ′ dx" dy" dx dy D x,2 (t) = t 0 ∂ x H t-s (x -x ′ , y -y)ρ (ε ′ ) (x ′ , y ′ , s)ρ (ε ′ ) (x", y", s) × δ (ε ′ ) (y ′ -y") sgn (ε) (x ′ -x") -sgn (ε ′ ) (x ′ -x") dx ′ dy ′ dx" dy" dx dy D y,1 (t) = t 0 ∂ x H t-s (x -x ′ , y -y ′ )ρ (ε ′ ) (x ′ , y ′ , s)ρ (ε ′ ) (x", y", s) × sgn (ε) (y ′ -y") δ (ε) (x ′ -x") -δ (ε ′ ) (x ′ -x") dx ′ dy ′ dx" dy" dx dy D y,2 (t) = t 0 ∂ y H t-s (x -x ′ , y -y)ρ (ε ′ ) (x ′ , y ′ , s)ρ (ε ′ ) (x", y", s) × δ (ε ′ ) (x ′ -x") sgn (ε) (y ′ -y") -sgn (ε ′ ) (y ′ -y") dx ′ dy ′ dx" dy" dx dy.
In order to study D x,1 , we make use of the following quantity

∂ x H t-s (x -x ′ , y -y ′ )ρ (ε ′ ) (x ′ , y ′ , s)I ε,ε ′ (x ′ , y ′ , s) dx ′ dy ′ with I ε,ε ′ (x ′ , y ′ , s) = ρ (ε ′ ) (x", y", s)sgn (ε) (x ′ -x") δ (ε) (y ′ -y") -δ (ε ′ ) (y ′ -y") dx" dy". Since δ (ε) (u) = 1 2 d
du sgn (ε) (u), the latter can be rewritten by integrating by parts

I ε,ε ′ (x ′ , y ′ , s) = 1 2 ∂ y ρ (ε ′ ) (x", y", s)sgn (ε) (x ′ -x") sgn (ε) (y ′ -y") -sgn (ε ′ ) (y ′ -y") dx" dy".
The Cauchy-Schwarz inequality yields

|I ε,ε ′ (x ′ , y ′ , s)| ≤ 1 2 (1 + |x"| 2 ) ∂ y ρ (ε ′ ) (x", y", s) 2 dx" dy" 1/2 × dx" 1 + |x"| 2 sgn (ε) (y ′ -y") -sgn (ε ′ ) (y ′ -y") 2 dy" 1/2 ≤ √ π 2 (1 + |z"| 2 )|∇ρ (ε ′ ) (z", s)| 2 dz" 1/2 ∆ ε,ε ′ (y ′ ), where ∆ ε,ε ′ (y ′ ) = sgn (ε) (y ′ -y") -sgn (ε ′ ) (y ′ -y") 2 dy" = 2 π y ′ -y" 0 e -v 2 2ε 2 dv ε - y ′ -y" 0 e -v 2 2|ε ′ | 2 dv ε ′ 2 dy" = 2 π u/ε u/ε ′ e -v 2 /2 dv 2 du.
In particular this quantity does not depend on y ′ . Clearly, for any fixed u ∈ R, we have lim

ε,ε ′ →0 u/ε u/ε ′ e -v 2 /2 dv = 0.
Furthermore, for 0 < ε, ε ′ ≪ 1, it can be dominated as follows

u/ε u/ε ′ e -v 2 /2 dv = u/ε u/ε ′ e -v 2 /4 e -v 2 /4
dv ≤ e -u 2 /2 e -v 2 /4 dv which lies in L 2 (R). Therefore the Lebesgue theorem tells us that lim

ε,ε ′ →0 ∆ ε,ε ′ = 0.
We go back to D x,1 that we split into

D x,1 (t) = t-η 0 ... ds + t t-η
... ds with 0 < η ≪ t ≤ T < ∞ to be determined. The integral on (0, t -η) can be estimated owing to the previous manipulations and the Cauchy-Schwarz inequality; we get

t-η 0 ... ds ≤ ρ (ε) L ∞ t-η 0 C 0 √ t -s sup x ′ ,y ′ |I ε,ε ′ (x ′ , y ′ , s)| ds ≤ C 0 √ π 2 ρ (ε) L ∞ ∆ ε,ε ′ t-η 0 ds t -s 1/2 t-η 0 (1 + |z| 2 )|∇ρ (ε) (z, s)| 2 ds dz 1/2 ≤ C T ∆ ε,ε ′ ln(t/η)
for a certain C T > 0, that comes from the estimates in Lemma 4.4. For the integral over (t -η, t), we claim that we can find a constant, still denoted

C T > 0, such that t t-η ... ds ≤ t t-η C 0 √ t -s δ (ε) (y ′ -y") + δ (ε ′ ) (y ′ -y") (38)
×ρ (ε) (x ′ , y ′ , s)ρ (ε ′ ) (x", y", s) dx ′ dy dx" dy" ds ≤ C T √ η.
This conclusion follows from uniform bounds (with respect to ε, ε ′ and s) of expressions like

J ε,ε ′ (s) = δ (ε) (y ′ -y")ρ (ε) (x ′ , y ′ )ρ (ε ′ ) (x", y") dz ′ dz". Let us set ρ(ε) (x", y ′ , s) = δ (ε) (y ′ -y")ρ (ε) (x", y", s) dy".
We control J ε,ε ′ (s) by using moments. Indeed, we get

J ε,ε ′ (s) = ρ (ε) (x ′ , y ′ , s)ρ (ε ′ ) (x", y ′ , s) dx ′ dy ′ dx" ≤ 1 2 1 + x ′2 1 + x" 2 |ρ (ε) | 2 (x ′ , y ′ , s) dx ′ dy ′ dx" + 1 2 1 + x" 2 1 + x ′2 |ρ (ε ′ ) | 2 (x", y ′ , s) dx ′ dy ′ dx" ≤ π 2 (1 + x 2 )|ρ (ε) | 2 (x, y, s) dx dy + π 2 (1 + x 2 )|ρ (ε ′ ) | 2 (x, y, s) dx dy.
Owing to Lemma 4.4, we already know that the first integral in the right hand side is bounded (the constant depends on the final time). For the second term, we simply write

(1 + x 2 )|ρ (ε ′ ) | 2 (x, y, s) dx dy ≤ (1 + x 2 ) δ (ε) (y -y ′ ) δ (ε) (y -y ′ )ρ (ε) (x, y ′ , s) dy ′ 2 (x, y, s) dx dy ≤ (1 + x 2 ) δ (ε) (y -y ′ ) dy ′ × δ (ε) (y -y ′ )|ρ (ε) (x, y ′ , s)| 2 dy ′ dx dy ≤ (1 + x 2 )|ρ (ε) (x, y ′ , s)| 2 δ (ε) (y -y ′ ) dy dx dy ′ ≤ (1 + x 2 )|ρ (ε) (x, y ′ , s)| 2 dx dy ′
which is thus also bounded uniformly with respect to ε, ε ′ > 0 and 0 ≤ s ≤ T < ∞. Finally, we arrive at

|D x,1 (t)| ≤ C T ln(t/η) ∆ ε,ε ′ + √ η which holds for any 0 < η ≪ t ≤ T < ∞. It shows that lim (ε,ε ′ )→0 D x,1 (t) = 0 uniformly on [0, T ].
The analysis of D x,2 is simpler; it relies on the following observation

δ (ε ′ ) (y ′ -y") sgn (ε) (x ′ -x") -sgn (ε ′ ) (x ′ -x") ρ (ε ′ ) (x", y", s) dx" dy" ≤ ρ (ε ′ ) ∞ |sgn (ε) (x ′ -x") -sgn (ε ′ ) (x ′ -x")| dx" ≤ 2 π u/ε u/ε ′ e -v 2 /2 dv du = ∆ε,ε ′ .
A straightforward adaptation of the argument used for studying ∆ ε,ε ′ shows that lim ε,ε ′ →0 ∆ε,ε ′ = 0 and we have

|D x,2 (t)| ≤ t 0 C 0 √ t -s ρ (ε) (s, •) L 1 ∆ε,ε ′ ds ≤ C T ∆ε,ε ′
for any 0 ≤ t ≤ T < ∞. Of course, D y,1 and D y,2 can be dealt with in a similar manner.

Coming back to (32), we arrive at

(ρ (ε) -ρ (ε ′ ) )(t, •) L 1 ≤ A ε,ε ′ + A(R, λ) + B(R) t 0 (ρ (ε) -ρ (ε ′ ) )(s, •) L 1 √ t -s ds, (39)
which holds for any 0

≤ t ≤ T < ∞ and 0 < R < ∞ with A ε,ε ′ = sup 0≤t≤T A ε,ε ′ (t) + sup 0≤t≤T D ε,ε ′ (t), A(R, λ) = (β 2 + γ 2 ) 1 + λ λ 2 √ T e ϕ(λ)T e -αλR , B(R) = (β 1 + γ 1 )RM, (40)
with M = sup ε>0 ρ (ε) L ∞ , which is known to be finite. We should bear in mind the fact that β 2 and γ 2 depend on λ too, through the exponential moments E 0 (2λ(1 + α)). Applying the singular Grönwall Lemma 3.4 leads to

(ρ (ε) -ρ (ε ′ ) )(t, •) L 1 ≤ A ε,ε ′ + A(R, λ) E 1/2 B(R) 2 √ t .
We remind the reader that the Mittag-Leffler function is explicitely known

E 1/2 (z) = ∞ k=0 z k Γ(1 + k/2) = e z 2 erfc(-z) = 2 √ π e z 2 z
-∞ e -u 2 du.

We are paying attention to the term

A(R, λ)E 1/2 B(R) 2 √
t . This is where we make use of (31)

to control E 0 (2λ(1 + α)) in the coefficients β 2 , γ 2 . As far as λ ≥ 1, we have ϕ(λ) ≤ 4D(1 + α) 2 λ 2 . Therefore, up to some irrelevant constant hereafter denoted by K > 0, the quantity of interest can be dominated by

√ T λ exp (DT + p 2 )4(1 + α) 2 λ 2 -αRλ + p 0 + (β 1 + γ 1 ) 2 M 2 4 T R 2 .
The exponent recasts as

4(DT +p 2 )(1+α) 2 λ- αR 8(DT + p 2 )(1 + α) 2 2 -R 2 α 2 16(DT + p 2 )(1 + α) 2 - (β 1 + γ 1 ) 2 M 2 4 T +p 0 .
We start by picking 0 < T < T ⋆ small enough, so that

α 2 16(DT + p 2 )(1 + α) 2 - (β 1 + γ 1 ) 2 M 2 4 T ≥ α 2 8p 2 (1 + α) 2 = q 2 > 0 holds for any 0 ≤ t ≤ T ⋆ . Next, let ω > 0. We can find R = R(ω) large enough so that K √ T e p0 e -R 2 q2 ≤ ω 2
holds. Possibly enlarging R(ω), we also suppose that

αR 8(DT + p 2 )(1 + α) 2 ≥ 1.
We then make use of the estimates with

λ = αR 8(DT + p 2 )(1 + α) 2 which leads to A(R, λ)E 1/2 B(R) 2 √ t ≤ ω 2 .
Finally, there exists ε(ω) > 0 small enough such that for any 0 < ε, ε ′ ≤ ε(ω) we get

A ε,ε ′ E 1/2 B(R) 2 √ t ≤ ω 2 . It follows that (ρ (ε) -ρ (ε ′ ) )(t, •) L 1 ≤ ω holds for any 0 ≤ t ≤ T ≤ T ⋆ , provided 0 < ε, ε ′ ≤ ε(ω).
We extend this result on any time interval by repeating the reasoning on subintervals of length smaller than T ⋆ . Therefore ρ (ε) ε>0 is a Cauchy sequence in the Banach space C([0, T ], L 1 (R 2 )) and it converges strongly to a solution of ( 1)-( 3). The proof can be readily adapted to establish the uniqueness of the solution of ( 1)-( 3) for a symmetric initial data verifying (31).

A convergence rate for ρ (ε)

Following the same strategy as in the proof of Theorem 4.5, it is possible to give a rate of convergence for ρ (ε) .

Theorem 4.10 Let T be a fixed time. Let ρ (ε) be the symmetric solution of (5)-( 6) with initial data ρ 0 (ρ 0 is assumed to be symmetric), and let ρ be the symmetric solution of (1)-( 2) with same initial data ρ 0 . Then there exist constants C(ρ 0 , T ) and 0 < ν(ρ 0 , T ) < 1 depending on both ρ 0 and T , such that sup

t∈[0,T ] ||(ρ (ε) -ρ)(t)|| L 1 ≤ C(ρ 0 , T )ε 1 2 ν(ρ0,T ) (41)
Remark 4.11 Observe that ν(ρ 0 , T ) is always smaller than 1, and it has the following asymptotic behavior lim

T →0 ν(ρ 0 , T ) = 1, lim T →+∞ ν(ρ 0 , T ) = 0,
for any ρ 0 . Note the 1/2 factor: with the present proof the convergence rate cannot be better than

ε 1/2 .
Proof. The idea is to revisit the computations in Section 4.2.2, in order to estimate more accurately the distance between ρ (ε) and ρ, solution of the singular PDE. Since we have used estimates that are uniform with respect to ε, we may simply take ε ′ = 0 in the computations performed above. It leads to the following observations:

• A term: We take the same initial condition for ρ (ε) and ρ, hence the error related to the initial condition simply vanishes: A ε (t) = 0.

• B and C terms: We use Lemmas 4.7 and 4.8, with ε ′ = 0.

• D term: We need to estimate

∆ ε = 4 π ∞ 0 ∞ u/ε e -v 2 /2 dv 2 du. Since v 2 /2 ≥ x 2 /2 + x(v -x), we get ∞ x e -v 2 /2 dv ≤ e -x 2 /2 ∞ 0 e -xs ds ≤ e -x 2 /2 x .
Thus, for any α > 0 we obtain The estimate for D x,1 (t) reads, for any 0 < η < t,

∆ ε ≤ 4 π α 0 π 2 du + ∞ α ε 2 u 2 e -u 2 /ε 2 du ≤ 4 π α π 2 + ε ∞ α/ε 1 s 2 e -
|D x,1 (t)| ≤ C T ∆ ε ln(t/η) + √ η .
Choosing η = ∆ ε (it is possible to do marginally better), we obtain, at the price of modifying

C T , |D x,1 (t)| ≤ C T ∆ ε ln(t/∆ ε ) ,
which, according to the above estimate for ∆ ε , yields

|D x,1 (t)| ≤ C T √ ε ln(t/ε).
Since D x,2 ≤ C T ε, we see that D x,1 is the largest contribution to D ε .

We use now (39)-(40) of the previous section with ε ′ = 0:

||(ρ (ε) -ρ)(t)|| L 1 ≤ (A ε + Ã(R, λ))E 1/2 (β 1 + γ 1 ) 2 M 2 R 2 T 4 . (42)
The contribution to A ε coming from the initial condition vanishes, since we choose the same initial condition for ρ (ε) and ρ. The second contribution to A ε comes from the "D terms", which are smaller than C T √ ε ln(T /ε). We can play the same game as in the proof of the Cauchy property: write E 1/2 (z) ≤ ce z 2 for some c, and observe that the exponent in (42) can be rewritten as

4(DT +p 2 )(1+α) 2 λ - αR 8(DT + p 2 (1 + α) 2 2 -R 2 α 2 4(DT + p 2 )(1 + α) 2 - (β 1 + γ 1 ) 2 M 2 T 4 +p 0 .
(43) We choose T = T * small enough so that the second term, proportional to R 2 is negative, which means

α 2 4(DT * + p 2 )(1 + α) 2 - (β 1 + γ 1 ) 2 M 2 T * 4 = q 2 (T * ) > 0
Then we choose λ such that the first term in (43) vanishes. We finally obtain sup

t∈[0,T * ] ||(ρ (ε) -ρ)(t)|| L 1 ≤ CA ε exp (K T * R 2 ) + C ′ exp (-q 2 R 2 ) ( 44
)
where C and C ′ depend on T * , and

K T = (β 1 + γ 1 ) 2 M 2 T 4 .
We now choose R to minimize the right hand side of (44). For instance, taking R such that exp [-(K

T * + q 2 (T * ))R 2 ] = A ε yields, for a modified C, sup t∈[0,T * ] ||(ρ (ε) -ρ)(t)|| L 1 ≤ C √ ε ln(T * /ε) ν . ( 45
) with ν = q 2 (T * ) K T * + q 2 (T * )
Slightly decreasing ν to absorb the logarithmic term, this proves the claim for any T < T * . For T > T * , we divide [0, T ] into subintervals of size T * , and apply the previous strategy for each subinterval. We have to take into account the error related to initial condition at the beginning of each subinterval. This error is given by the total error at the end of the previous subinterval. Thus we have to reintroduce an error related to initial data. Calling E k the bound on the error at the end of the interval [(k -1)T * , kT * ], and A (k) ε the A ε term to be considered on the interval [kT * , (k + 1)T * ], we have

A (k) ε ≤ C T * √ ε ln(T * /ε) + C T * E k ≤ C T * E k ,
where C T * can take different values, but remains a constant depending on ρ 0 , T * , and not on ε.

With the same reasoning as above, we conclude with

E k+1 ≤ C T * E ν k .
Since T * is of order 1, we have to repeat the argument on a finite number of subintervals to reach the prescribed time T . Each iteration of course decreases the convergence rate, and increases the prefactor, but for any T , we can guarantee a finite ν, as claimed.

Particle approximation

We consider now an N -particle description of the dynamics. Namely, let Z

(ε) i = (X (ε) i , Y (ε)
i ) be the solution of the stochastic differential system dX

(ε) i,t = 1 N j =i K (ε) x (Z (ε) i,t -Z (ε) j,t ) dt + √ 2D dB i,x,t , (46) dY
(ε) i,t = 1 N j =i K (ε) y (Z (ε) i,t -Z (ε) j,t ) dt + √ 2D dB i,y,t , (47)
where B i,x and B i,y are independent Brownian motions. Here and below, the interaction kernel is given by

K (ε) x (z) = -sgn (ε) (x)δ (ε) (y) K (ε) y (z) = -sgn (ε) (y)δ (ε) (x),
with z = (x, y). It is then clear that K (ε) x Lip = C/ε 2 , and the same holds true for K (ε) y . We assume that the initial conditions for the particles' trajectories

Z (ε) i,t t=0 = Z (ε) i,0
are independent random variables, with common law ρ 0 . In the discussion, we naturally assume that ρ 0 is a probability density. Accordingly, for both ρ and ρ (ε) solutions of (1) and ( 5) respectively, associated to the initial data ρ 0 , we have

ρ dz = ρ (ε) dz = ρ 0 dz = 1.
Moreover, we assume throughout this section that ρ 0 is such that the symmetric existence theorem works as we shall use the rate of convergence established in this framework. We associate to the solutions of this system (46)-(47), the empirical measure

ρ (ε),N = 1 N N i=1 δ(z -Z (ε) i ).
Note that the interaction force in ( 46)-(47) has been rescaled by the 1/N factor (roughly speaking we have replaced the kernel K (ε) by 1 N K (ε) ), so that the total force exerted on a given particle remains of order 1; this is the so-called mean field regime. We refer the reader to the surveys [START_REF] Bolley | Limite de champ moyen de systèmes de particules[END_REF][START_REF] Golse | The mean-field limit for the dynamics of large particle systems[END_REF] for an introduction to such regimes. The goal of this section is to investigate the convergence of this particle approximation to ρ, the solution of the singular PDE (1) in the regime N → ∞, ε → 0.

The analysis uses the Wasserstein distance, see [START_REF] Dobrušin | Vlasov equations[END_REF][START_REF] Villani | Optimal transport, old and new[END_REF] for a thorough discussion on this notion. The Wasserstein distance W 1 (µ, ν) between two probability measures µ, ν on R 2 is defined as

W 1 (µ, ν) = sup ϕ(z)µ(dz) -ϕ(z)ν(dz) , ϕ Lip ≤ 1 ,
where

ϕ Lip = sup x =y, x,y∈R 2 |ϕ(x) -ϕ(y)| |x -y| .
Note that W 1 determines the topology of tight convergence on the space of probability measures on R 2 , see [START_REF] Villani | Optimal transport, old and new[END_REF]Chap. 6].

Wasserstein metric is well defined on the set of probability measures with finite first moment. This is the case for ρ, the solution of the original PDE (1), see Proposition 2.1 as well as ρ (ε) , the solution of the regularized PDE [START_REF] Brézis | Analyse fonctionnelle. Théorie et applications[END_REF], see Proposition 2.2. It also holds true for the particle approximations ρ (ε),N (they are finite sums of Dirac delta distributions).

It turns out that W 1 is a well adapted tool to investigate the limit N → ∞, see [START_REF] Bolley | Limite de champ moyen de systèmes de particules[END_REF][START_REF] Dobrušin | Vlasov equations[END_REF][START_REF] Golse | The mean-field limit for the dynamics of large particle systems[END_REF][START_REF] Sznitman | Topics in propagation of chaos[END_REF]]. The strategy is to write

W 1 ( ρ (ε),N , ρ) ≤ W 1 ( ρ (ε),N , ρ (ε) ) + W 1 (ρ (ε) , ρ),
where ρ (ε) is the solution of the regularized PDE [START_REF] Brézis | Analyse fonctionnelle. Théorie et applications[END_REF]. The second term is controlled by the rate of convergence established in the previous section, and the first one by adapting "standard" MacKean-Vlasov estimates, as we are going to detail now. According to [START_REF] Sznitman | Topics in propagation of chaos[END_REF], we start by introducing an auxiliary system of interacting particles. The solution ρ (ε) of the regularized PDE is also the law of the solution of the system of SDE d X(ε

) i = (K (ε) x ⋆ ρ (ε) )( Z(ε) i ) dt + √ 2D dB i,x , (48) d
Ỹ (ε) i = (K (ε) y ⋆ ρ (ε) )( Z(ε) i ) dt + √ 2D dB i,y . (49)
Note that both Z

(ε) i = (X (ε) i , Y (ε)
i ) and Z(ε

) i = ( X(ε) i , Ỹ (ε)
i ) are driven by the same Brownian motions and we choose them to have the same initial condition. The system of stochastic differential equations (48)-(49) (respectively (46)-( 47)) has a unique solution, as the coefficients (t, z) → K

(ε) x ⋆ ρ (ε) (z) and (t, z) → K (ε)
y ⋆ ρ (ε) (z) are Lipschitz with respect to z and continuous with respect to t. Moreover, the law µ

(ε) of Z(ε) i = ( X(ε) i , Ỹ (ε) i ) is a (weak) solution of ∂ t µ (ε) = ∇ • (-K (ε) ⋆ ρ (ε) )µ (ε) + D∆µ (ε) µ (ε) t=0 = ρ 0 .
Since this equation has a unique solution, and ρ (ε) is a solution, it follows that µ (ε) = ρ (ε) . We define ρ (ε),N to be the empirical measure associated with the Z(ε) i :

ρ (ε),N = 1 N N i=1 δ(z - Z(ε) i ).
The following statement is an immediate corollary of Theorem 1 in [START_REF] Fournier | On the rate of convergence in Wasserstein distance of the empirical measure[END_REF]:

Proposition 5.1 Let 0 < T < ∞. Assume that there exist q > 2 and a constant C = C(T ), that depends on T but is independent of ε, such that

sup t∈[0,T ] |z| q ρ (ε) (dz) ≤ C. (50)
Then there exists a constant C = C(T ) independent of ε such that

sup t∈[0,T ] E[W 1 ( ρ (ε),N , ρ (ε) )] ≤ C √ N log (1 + N ). (51)
The uniform bound (50) holds true if, for example, the initial measure ρ 0 has a finite third moment. To prove this, one uses an argument similar to that in Proposition 2.2-iv (in effect one uses the same proof as that used for the a priori bound deduced for the original measure ρ in Proposition 2.1-iv). We state and prove now the main result of the section. Theorem 5.2 Let 0 < T < ∞ be a fixed time. Under the same conditions as in Theorem 4.10, we have

sup t∈[0,T ] E[W 1 ( ρ (ε),N , ρ)] ≤ Ce 2CT ε 2 √ N log (1 + N ) + C ρ ε 1 2 νρ , (52)
where C = C(T ) is the constant defined in Proposition 5.1, C is the Lipschitz constant of ε 2 K (ε) and C ρ = C(ρ 0 , T ), respectively, ν ρ = ν(ρ 0 , T ) are the constants arising from Theorem 4.10.

In particular, for any δ ∈ [0, 1 2 ), there exists ε = ε(δ, N ) and a constant Cρ = Cρ (δ) independent of N such that sup

t∈[0,T ] E[W 1 ( ρ (ε),N , ρ)] ≤ Cρ (log(N )) -1 4 νρ (53)
for any N ≥ 1.

Proof. Following Theorem 4.10, to establish (52) it suffices to prove that sup

t∈[0,T ] E[W 1 ( ρ (ε),N , ρ (ε) )] ≤ Ce 2CT ε 2 √ N log (1 + N ). ( 54
) Since both Z (ε) i = (X (ε) i , Y (ε)
i ) and Z(ε

) i = ( X(ε) i , Ỹ (ε)
i ) are driven by the same Brownian motions and have the same initial condition, we have

d dt (Z (ε) i,t - Z(ε) i,t ) = 1 N j K (ε) (Z (ε) i,t -Z (ε) j,t ) -(K (ε) ⋆ ρ (ε) )( Z(ε) i,t ) = (K (ε) ⋆ ρ (ε),N )(Z (ε) i,t ) -(K (ε) ⋆ ρ ε )( Z(ε) i,t ) = [K (ε) ⋆ ( ρ (ε),N -ρ (ε) )(Z (ε) i,t )] + [(K (ε) ⋆ ρ (ε) )(Z (ε) i,t ) -(K (ε) ⋆ ρ (ε) )( Z(ε) i,t )].
We note that K

(ε) (• -Z (ε) i,t ) is a function with Lipschitz constant less than C ε 2 . Hence K (ε) ⋆ ( ρ (ε),N -ρ (ε) )(Z (ε) i,t ) ≤ C ε 2 W 1 ( ρ (ε),N , ρ (ε) ).
Furthermore, using that K (ε) (z -•) is C/ε 2 -Lipschitz, and ρ ε dz = 1, we get

K (ε) ⋆ ρ (ε) (Z (ε) i,t ) -K (ε) ⋆ ρ (ε) ( Z(ε) i,t ) ≤ C ε 2 |Z (ε) i,t - Z(ε) i,t |. ( 55
) Then d dt |Z (ε) i,t - Z(ε) i,t | ≤ C ε 2 |Z (ε) i,t - Z(ε) i,t | + C ε 2 W 1 ( ρ (ε),N , ρ (ε) ).
Hence, since Z (ε) and Z(ε) share the same initial data, we arrive at

e -Ct ε 2 |Z (ε) i,t - Z(ε) i,t | ≤ t 0 Ce -Cs ε 2
ε 2 W 1 ( ρ (ε),N , ρ (ε) )(s) ds.

(56)

Now we write e -Ct ε 2 W 1 ( ρ (ε),N , ρ (ε) ) ≤ e -Ct ε 2 W 1 ( ρ (ε),N , ρ (ε),N ) + e -Ct ε 2 W 1 ( ρ (ε),N , ρ (ε) )

≤ 1 N N i=1 e -Ct ε 2 |Z (ε) i,t - Z(ε) i,t | + e -Ct ε 2 W 1 ( ρ (ε),N , ρ (ε) ) ≤ t 0 Ce -Cs ε 2 ε 2 W 1 ( ρ (ε),N , ρ (ε) )(s) ds + e -Ct ε 2 W 1 ( ρ (ε),N , ρ (ε) ),
where we have used first the triangle inequality, then a direct inequality for the W 1 distance between the two empirical measures, and finally (56). By taking the expectation and using (51), we obtain e -Ct ε 2 EW 1 ( ρ (ε),N , ρ which gives (53).

Numerical illustrations

The goal of this section is two-fold:

1. Illustrate the existence Theorems 4.1 and 4.5, and show that the minimal value for the diffusion we have identified in the statement is not optimal: the solution can apparently be global in time for D < 2C 2 M 0 ;

2. Illustrate the convergence for the particles approximation, and show that the actual rate of convergence as a function of N seems to be much better than suggested by Theorem 5.2.

For this purpose, we use a finite volume method introduced in [START_REF] Carrillo | A finite-volume method for nonlinear nonlocal equations with a gradient flow structure[END_REF] to study drift-diffusion equations with gradient structure. Of course, there is no gradient structure in the present case, but the method can be adapted and it is proved to be robust. Let us briefly explain the principles of the approach. We work on a Cartesian grid, with space steps ∆x, ∆y > 0. Given the time step ∆t > 0, we wish to update the numerical unknown with a finite volume formula which looks like ρ n+1 i,j = ρ n+1 i,j -∆t ∆x (F i+1/2,j -F i-1/2,j ) -∆t ∆y (G i,j+1/2 -G i,j-1/2 )

where we need to find a relevant definition for the numerical fluxes F, G. To this end, we rewrite the right hand side of (1) as

∇ • ρ(∇ ln(ρ) -F [ρ]) = ∂ x ρ(∂ x ln(ρ) + ∂ x U ) + ∂ y ρ(∂ y ln(ρ) + ∂ y V )
where U, V are the scalar functions defined by U (x, y, t) = |x -x ′ |ρ(x ′ , y, t) dx ′ , V (x, y, t) = |y -y ′ |ρ(x, y ′ , t) dx ′ .

We shall therefore apply the ideas in [START_REF] Carrillo | A finite-volume method for nonlinear nonlocal equations with a gradient flow structure[END_REF] directionwise. The flux F i+1/2,j is given by applying the upwinding principle with the "velocity" ξ = ∂ x ln(ρ) + ∂ x U which leads to F i+1/2,j = ξ i+1/2,j + ρ i,j + ξ i+1/2,j + ρ i+1,j .

The interface value is obtained by the mere centered difference ξ i+1/2,j = 1 ∆x ln(ρ i+1,j ) -ln(ρ i,j ) + U i+1,j -U i,j ,

where the integral that defines U can be evaluated by a quadrature rule (the rectangle rule, say). A similar construction applies to construct the flux G. The accuracy of the method can be improved by using a polynomial reconstruction of the density, with a suitable slope limiter, instead of the mere upwind scheme, in the spirit of the design of MUSCL schemes. We refer the reader to [START_REF] Carrillo | A finite-volume method for nonlinear nonlocal equations with a gradient flow structure[END_REF] for further details and the analysis of this scheme for gradient-flow equations. We can equally use a second-order Runge-Kutta method for the time integration. We do not explicitly introduce a regularization for the singular forces (2) in the code; we simply compute (2) by summing over rows or columns of the square grid. This corresponds to an effective regularization of the order of the grid spacing (typically ∆x = ∆y = 0.05 in the simulations presented below). For the particles simulations, we integrate directly the regularized equations ( 46)-(47) by using the Euler method. We typically use ǫ = 0.1. Fig. 1 shows a contour plot of ρ at late times for D = 0.15 obtained by using the finite volume method introduced in [START_REF] Carrillo | A finite-volume method for nonlinear nonlocal equations with a gradient flow structure[END_REF] (left plot) and the (mollified) particles approximation (right plot). Fig. 2 shows the evolution of the L 2 and L ∞ norms for various values of D. D = 0.15 is smaller than 2C 2 , the threshold of Theorem 4.1 (here M 0 = 1): the L 2 norm is not monotonically decreasing, but there is apparently no finite time singularity. Fig. 3 shows that particles simulations are reasonably close to the PDE simulations already for a number of particles much smaller than that suggested by Theorem 5.2.

Lemma 4 . 4 (

 44 Weighted L 2 estimates) Let ρ (ε) 0 ε>0 be a sequence of non negative functions bounded in L 1 ∩ L ∞ (R 2 ). Assume that ρ (ε) 0 satisfies (28) and sup ε>0 |z| 4 ρ (ε)

  s 2 ds . Choosing α = ε, this relation yields ∆ ε ≤ Cε where C is an absolute constant. A very similar reasoning applied to ∆ε = e -v 2 /2 dv du yields ∆ε ≤ Cε where again C is an absolute constant.

  (ε) )(t) ≤ (1 -e -Ct ε 2 ) C √ N log (1 + N ) + C ε 2 t 0 e -Cs ε 2 EW 1 ( ρ (ε),N , ρ (ε) )(s) ds. Cs ε 2 EW 1 ( ρ (ε),N , ρ (ε) )(s) ds.By the standard Grönwall's lemma we deduce thate -Ct ε 2 EW 1 ( ρ (ε),N , ρ (ε) )(t) 54). Using the triangle inequality and Theorem 4.10, (54) leads to (52). Moreover observe that for δ ∈ [0,1 2 ) and ε = ( 1-2δ 4CT log(N ))

Figure 1 :

 1 Figure1: Contour plot of the density ρ for D = 0.15, at time t = 5. The left plot is done using the finite volume method introduced in[START_REF] Carrillo | A finite-volume method for nonlinear nonlocal equations with a gradient flow structure[END_REF]. The right plot is done using the (mollified) particles approximation with 10 samples of 10 4 particles. Note that the noise due to the finite number of particles is still visible.

Figure 2 :

 2 Figure 2: L 2 norm (left) and L ∞ norm (right) as a function of time for D = 0.15 (black), D = 0.25 (red) and D = 0.35 (blue).

  We need to strengthen this local property to a global statement. The moment estimate and the L ∞ estimate in Lemma 3.3-i) and ii) respectively, allow us to justify that

	T
	0	|z|≥R

we can appeal to the Aubin-Simon lemma, see

[START_REF] Simon | Compact sets in L p (0, T ; B). Ann. Mat. Pura e Applicada[END_REF] Cor. 4, Sect. 8]

, to deduce that ρ n n∈N is relatively compact in L 2 ((0, T ) × B(0, R)) for any 0 < R < ∞.

The constant Cp is the constant appearing in the Gagliardo-Nirenberg-Sobolev inequality[START_REF] Dalibard | Laser cooling of an optically thick gas: the simplest radiation pressure trap ?[END_REF].

Acknowledgements

We are gratefully indebted to Nicolas Fournier and José Antonio Carrillo for many motivating discussions and helpful advices. J. Barré acknowledges the support of CNRS and Imperial College London which has made a visit of several months at the Math Department of ICL possible.


A Drogoul
email: audric.drogoul@inria.fr

G Aubert

B Cessac

P Kornprobst

A new nonconvex variational approach for sensory neurons receptive field estimation

Determining the receptive field of a visual sensory neuron is crucial to characterize the region of the visual field and the stimuli this neuron is sensitive to. We propose a new method to estimate receptive fields by a nonconvex variational approach, thus relaxing the simplifying and unrealistic assumption of convexity made by standard approaches. The method consists of solving a relaxed discrete energy minimization problem using a proximal alternating algorithm. We compare our approach with the classical spike-triggered-average technique on simulated data, considering a typical retinal ganglion cell as ground truth. Results show a high improvement in terms of accuracy and convergence with respect to the duration of the experiment.

Introduction

This paper deals with the estimation of receptive fields of individual visual neurons. Knowing the receptive field of a particular neuron is a crucial information to understand which region of the visual field and which stimulus the neuron is sensitive to. Thus, biologists always devote efforts in their experimental protocol to characterize the neurons they are recording.

The goal is to characterize the relation between the stimulus and the neuron response which is a set of action potentials (also called spikes). To do so, it is classical to assume that response follows a linear-nonlinear Poissonian (LNP) model [START_REF] Chichilnisky | A simple white noise analysis of neuronal light responses[END_REF][START_REF] Mcfarland | Inferring Nonlinear Neuronal Computation Based on Physiologically Plausible Inputs[END_REF][START_REF] Paninski | Maximum likelihood estimation of cascade point-process neural encoding models[END_REF] (Fig. 1): given a visual stimulus s(x, t) : Ω ⊂ R 2 × [0, T ] → R, where Ω is the spatial domain and T is the duration of the experiment, the neuron generates a sequence of n(T ) spikes times {t i } 1≤i≤n(T ) such that {t i } 1≤i≤n(T ) is generated by a Poisson process of rate r(t) = S ((s × u)(t)) , (P)

where S(.) is a nonlinear function and u : Ω d = Ω × [-d, 0] -→ R is the so-called receptive field which corresponds to the linear part of the processing with d > 0 the length of its temporal support.

The operator × is defined by (s × u)(t) = 0 -∞ s(., t + τ ), u(., τ ) Ω dτ, where < ., . > Ω denotes the inner product on Ω. Given hypothesis (P) the problem of estimating the receptive field u(x, t) is an inverse problem: given a stimulus s(x, t) of duration T and the n(T ) spikes {t i } 1≤i≤n(T ) from the sensory neuron responding to s(x, t), recover the unknown receptive field u(x, t).

This inverse problem can be formulated using a Bayesian approach. Using the Bayes rule, a classical method is to search for the receptive field u(x, t) as minimum of the anti-log of the a posteriori probability (see, e.g., [START_REF] Paninski | Maximum likelihood estimation of cascade point-process neural encoding models[END_REF][START_REF] Pillow | Likelihood-based approaches to modeling the neural code. Bayesian brain: Probabilistic approaches to neural coding[END_REF][START_REF] Park | Receptive Field Inference with Localized Priors[END_REF]): Here the stimulus is a white noise sequence which is the classical stimulus used for the spike-triggered averaged approach (STA, see [START_REF] Chichilnisky | A simple white noise analysis of neuronal light responses[END_REF]).

inf u -log ρ(u|{t i } 1≤i≤n(T ) ) = inf u -log(ρ({t i } 1≤i≤n(T ) |u)) -log(ρ(u)) , (1)
where ρ(u) denotes the probability density of the random variable u. Then, assuming that ρ(u) is a distribution of the form ρ(u) ∝ e -J(u) , one can show that (1) can be rewritten as [START_REF] Drogoul | A new nonconvex variational approach for sensory neurons receptive field estimation[END_REF]:

inf u E(u) = ψ(s × u) + J(u), (2)
where ψ(s × u) = -log(ρ({t i } 1≤i≤n(T ) |u)) is the data fidelity term defined by

ψ(z) = T 0 S(z(τ ))dτ - n(T ) i=1 log(S(z(t i ))), (3)
with z(t) = (s × u](t), and J(u) is the prior term to infer qualitative properties to the solution, according to what is known about the general shape of a receptive field. To study this problem, it is classical to assume that the nonlinearity S(.) is a ramp function or an exponential. The computation of ( 2) is then simple since it is a convex problem to solve [START_REF] Paninski | Maximum likelihood estimation of cascade point-process neural encoding models[END_REF][START_REF] Pillow | Likelihood-based approaches to modeling the neural code. Bayesian brain: Probabilistic approaches to neural coding[END_REF][START_REF] Mcfarland | Inferring Nonlinear Neuronal Computation Based on Physiologically Plausible Inputs[END_REF]. On the opposite, in this paper, we consider the case when the nonlinearity S(.) is a sigmoid, which is more realistic from a physiological point of view; in particular it models the fact that the neuron has a bounded firing rate. However, since S(.) is nonconvex, the data fidelity term ψ(s × u) also becomes nonconvex. This is known to be harder to minimize numerically.

Concerning J(.), we choose it to impose two properties on the solution. The first is that u should be localized in space and time since neurons are sensitive to a particular region of the visual field. This can be imposed by a sparsity constraint term. Here we choose to use a convex relaxation of the sparsity and we penalise the L 1 -norm of u. The second is that u should be smooth. We propose that u should belong to the space BV 2 that contains piecewise linear functions. This kind of regularity constraint has also been previously used in the context of image restoration [START_REF] Bergounioux | A Second-Order Model for Image Denoising[END_REF]. This allows to recover functions with fast smooth variations more precisely than using a simple ∇u 2 L 2 -regularity constraint. So the prior term will be

J(u) = λ u L 1 (Ω d ) + µ|u| BV 2 (Ω d ) , (4)
where λ, µ > 0 are weights associated respectively to the sparsity and to the regularity of u.

So this paper is about the study of the problem ( 2) with ( 3)-( 4) In Sect. 2, we study the discrete version of the problem (2). We introduce a relaxed problem so that the solution can be computed more easily. We justify it theoretically and introduce an alternating minimizing algorithm converging toward a critical point of the relaxed energy. In Sect. 3, we test the approach on simulated data to provide a quantitative evaluation with comparisons to the classical spike triggered averaged technique (STA, see [START_REF] Chichilnisky | A simple white noise analysis of neuronal light responses[END_REF]). In Sect. 4, we describe future work.

Study of the discrete problem: Well-posedness and algorithm

Problem definition

In this section, we study a discrete version of (2). The stimulus is a sequence of N t images each presented during a period of ∆t so that the duration of the experiment is T = N t ∆t. Each image is of size N x × N y pixels. Receptive field is of size N x × N y × D where D is a fixed depth. We introduce the real vector spaces X , Y and Z so that:

u ∈ X = R Nx×Ny×D , s ∈ Y = R Nx×Ny×Nt , z ∈ Z = R Nt , (5)
endowed with the scalar product ., . and the associated norm . 2 = ., . . We denote by . 1 the l 1 -norm. We do the following hypothesis on the sigmoid function S(.).

Hypothesis 1. There exist c > 0 and θ 1 < θ 2 such that :

(i) S = 0 ] -θ 1 ], S increasing in [θ 1 , θ 2 ] and S = c on [θ 2 , +∞[; (ii) S ∈ C 1 (R) and analytic on [θ 1 , θ 2 ]; (iii) z → -log(S(z)) is convex on ]θ 1 , +∞[.
Given these notations, if we denote ξ = (ξ i ) 1≤i≤Nt the number of spikes per time step and if we assume that ∆t = 1 without loss of generality, then the data fidelity term (3) can be written as

ψ ξ (z) = Nt i=1 ψ ξ i (z i ), with ψ ξ i (z i ) = S(z i ) -ξ i log(S(z i )), (6)
where ψ ξ i (z) has a limited domain of definition denoted by

Z + ξ i equal to R if ξ i = 0, and ]θ 1 , +∞[ if ξ i > 0.
One can easily verify that the data fidelity term is nonconvex when ξ i ≤ c.

The prior term is approximated by

J ε (u) = λ u 1 + µI ε (u) with I ε (u) = (i,j,k) ε 2 + (Hu) i,j,k 2 2 , (7)
where the introduced parameter ε > 0, is small enough to make the prior term differentiable and H : X → X 9 is a matrix equal to the discrete Hessian. It is clear that by introducing a dual variable, by following for example [START_REF] Chambolle | [END_REF], we could have avoid to introduce this regularization.

We did not use this approach for computational time reason: we have to solve a primal dual problem at each iteration and compute a projection on the set {H p, p ∈ X 9 , |p i,j,k | ≤ 1} by a fixed point algorithm based on [START_REF] Chambolle | An Algorithm for Total Variation Minimization and Applications[END_REF]. Denoting by × d the discretization of × we get the discrete problems associated to (2):

inf u∈X E(u) = ψ ξ (s × d u) + i Z + ξ (s × d u) + λ u 1 + µI ε (u), (8)
where i Z + ξ denotes the characteristic function associated to Z + ξ = Nt i=1 Z + ξ i . In [START_REF] Kurdyka | On gradients of functions definable in o-minimal structures[END_REF], the data fidelity term depends only on s × d u and when minimizing w.r.t. u it is difficult to guarantee that s × d u remains in the domain of definition of ψ ξ . To solve this problem, we propose to introduce an auxiliary variable z and the following relaxed formulation:

inf z∈Z,u∈X E α (z, u) = ψ ξ (z) + i Z + ξ (z) + α 2 s × d u -z 2 2 + λ u 1 + µI ε (u), (9)
where the term in α penalizes the difference between z and s × d u. Another interest of this relaxed problem is that now the problem in z containing the nonconvexity is separable (it leads to N t one-dimensional independent problems) and we will show that the energy is separately convex w.r.t. to z and u under some assumptions.

Well-posedness

We start with the following two propositions which result from standard arguments:

Proposition 2.1 (Existence). Let α > 0 and λ > 0 (resp. λ > 0), the energy E α (z, u) (9) (resp. E(u) (8)) is coercive and admits at least a minimizer. The following proposition establishes the link between solutions of the relaxed problem and solutions of the initial one.

Proposition 2.3. Let α be a positive integer, then the sequence {(z α , u α )} α≥1 solution of the minimization problem (9) associated to E α (z, u) is bounded and (i) All its cluster points are couples (s × ū, ū) such that ū is a solution of (8).

(ii) The infimum converges:

inf (z,u) E α (z, u) = E α (z α , u α ) -→ α→+∞ inf u E(u). (iii) If E(u) admits a unique minimizer ū, then (z α , u α ) -→ α→+∞ (s × ū, ū).
Proof. The proof of this proposition uses the epi-convergence of E α (z, u) toward E(u) i {z=s× d u} (see [START_REF] Tyrrell Rockafellar | Variational Analysis, volume 317 of Grundlehren der mathematischen Wissenschaften[END_REF]).

Algorithm

To compute a solution of ( 9), we propose the proximal alternating minimization introduced in [START_REF] Attouch | Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the kurdyka-lojasiewicz inequality[END_REF] in a general context. Given an initial condition u 0 ∈ X (e.g., u 0 = 0), the algorithm consists of the following two steps:

z (k+1) ∈ arg min z E α (z, u (k) ) + 1 2β (k) z -z (k) 2 2 , (10a)
u (k+1) ∈ arg min u E α (z (k+1) , u) + 1 2γ (k) u -u (k) 2 2 , ( 10b
)
where β (k) , γ (k) are sequences of parameters belonging to [r -, r + ] with 0 < r -< r + for all k ≥ 0 (note that this is the only condition on these parameters to obtain convergence). Note that quadratic terms in [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF] are necessary to show the convergence and that the definition of ( 10) leads to a decrease of the energy verifying:

E α (z (k+1) , u (k+1) ) + 1 2β (k) z (k+1) -z (k) 2 2 + 1 2γ (k) u (k+1) -u (k) 2 2 ≤ E α (z (k) , u (k) ).
From a numerical point of view, this algorithm becomes easy to tackle:

• If we choose S(.) as a piecewise cubic function verifying (1), when ξ i = 0, then problem (10a) can be solved analytically thanks to its separability: it is equivalent to compute the proximal operator [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF] of S(.) up to a multiplicative constant depending on α and β (k) . When α + 1 β (k) >inf ess S (.) the problem is strictly convex so that the proximal operator is uni-valued. Otherwise, in the general case, we use a Newton algorithm.

• Problem (10b) is convex and we solve it by a standard forward-backward proximal algorithm of type Nesterov [START_REF] Beck | Fast Gradient-Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems[END_REF].

Theorem 2.1. (Convergence) Algorithm (10a)-(10b) generates a sequence (z (k) , u (k) ) k∈N which converges to a critical point (z, u) of E α (z, u) (i.e., 0 ∈ ∂E α (z, u)) and the sequence (z (k) , u (k) ) is l 1 (N).
Proof. Energy E α (z, u) verifies the Kurdyka-Lojasiewicz property (see [START_REF] Kurdyka | On gradients of functions definable in o-minimal structures[END_REF] and Definition 7 in [START_REF] Attouch | Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the kurdyka-lojasiewicz inequality[END_REF]) and is coercive from Proposition 2.1. By applying Theorem 9 of [START_REF] Attouch | Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the kurdyka-lojasiewicz inequality[END_REF] we get the result.

Numerical results

In this section we show quantitative results obtained using simulated data. We assume that we look for the receptive field of a neuron which can be approximated by a LNP model and that the nonlinearity is a known sigmoid (a piecewise cubic approximation verifying Hypothesis 1).

Here we choose as the ground truth a receptive field shape that corresponds to a class of retinal ganglion cells: the function u is separable in space and time, it is a difference of Gaussian in space and a difference of exponential in time. It is discretized as a spatio-temporal volume of size 20 × 20 × 30 pixels. This is illustrated in Fig. 2 (first row, sample temporal slices).

For the stimulus, we choose a sequence of uniform random binary images each presented during a time step ∆t = 1. Here we use images of size 20 × 20 pixels, with block size 4 × 4 pixels (as illustrated in Fig. 1). Using this stimulus allows to apply the STA technique [START_REF] Chichilnisky | A simple white noise analysis of neuronal light responses[END_REF] to make a comparison.

In Fig. 2 we show a comparison between STA [START_REF] Chichilnisky | A simple white noise analysis of neuronal light responses[END_REF] and our approach (β (k) = γ (k) = 10 ∀k, λ = 10, = 10 -8 , µ = 100, c = 33, θ 1 = -7, θ 2 = 107). This result is obtained using a stimulus of N t = 1000 images (giving around n(T ) = 500 spikes). Results show that the spatial resolution of STA is constrained by the size of the block of the stimulus, by definition of the STA (PSNR=22.1dB), while we obtain a much better estimation (PSNR=33.1dB). Concerning the convergence, different aspects are shown in Fig. 3. In Fig. 3(left), we show that our approach converges faster than STA as the number of spikes increases, which is directly related to the duration of the experiment. In other words, for a fixed duration of experiment, our approach provides a better estimate of the receptive field than STA. In Fig. 3(middle), the main observation is that error decays linearily (in log scale) with a slope of the convergence equal for both methods which can be explained by the central limit Theorem. Finally, Fig. 3(right) illustrates the speed of convergence of our algorithm.

GT

Conclusion

Overall, this study presents for the first time an approach to deal with the non-convex case allowing efficient extraction of receptive fields with great accuracy, as exemplified on synthetic data. Another advantage of this approach is that we are no longer constrained to use white noise as input as needed in STA, so that one can investigate what could be an optimal stimulus (e.g., [START_REF] Pamplona | A super-resolution approach for receptive fields estimation of neuronal ensembles[END_REF]) to further improve the results. Future work will focus on validating our approach on real cell recordings.

Figure 1 .

 1 Figure 1. Illustration of the LNP model (P). Two of S(.) are illustrated: ramp-or nonlinearities.Here the stimulus is a white noise sequence which is the classical stimulus used for the spike-triggered averaged approach (STA, see[START_REF] Chichilnisky | A simple white noise analysis of neuronal light responses[END_REF]).

Proposition 2 . 2 .

 22 If α >inf R ess S then problem E α (z, u) (9) is convex w.r.t to z and u separately (but not convex w.r.t. (z, u)).

Figure 2 .

 2 Figure 2. Example of result for simulated data.Each column shows a sample of temporal slices of (first line) the ground truth (GT), i.e., the receptive field to recover, (second line) the result obtained using the classical STA approach, and (third line) the result obtained using our approach. Last column shows temporal profile corresponding to a central position.

l 2 Figure 3 .

 23 Figure 3. Convergence of the approach: (left) covariance error (1 -cov(u,u 0 ) σuσu 0 ) as a function of the number of spikes in semi-log scale, (middle) l 2 -norm error ( uu 0 2 ) as a function of the number of spikes in log-log scale, and (right) l 2 -norm error as a function of the number of iterations for N t = 1000.

Acknowledgements This work was partially supported by the EC IP project FP7-ICT-2011-9 no. 600847 (RENVISION). We also thank Daniela Pamplona (Inria, Biovision team) as well as Evelyne Sernagor and Gerrit Hilgen (Institute of Neuroscience, Newcastle University, UK) for their helpful insights.