Alp Ár

Rich Árd

M Ész

Filippo Santambrogio

ADVECTION-DIFFUSION EQUATIONS WITH DENSITY CONSTRAINTS

Keywords: February 12, 2016 diffusive crowd motion model, Fokker-Planck equation, density constraints, optimal transportation. Mathematics Subject Classification: 35K61, 49D10, 49J45. 1

In the spirit of the macroscopic crowd motion models with hard congestion (i.e. a strong density constraint ρ ≤ 1) introduced by Maury et al. some years ago, we analyze a variant of the same models where diffusion of the agents is also taken into account. From the modeling point of view, this means that individuals try to follow a given spontaneous velocity, but are subject to a Brownian diffusion, and have to adapt to a density constraint which introduces a pressure term affecting the movement. From the PDE point of view, this corresponds to a modified Fokker-Planck equation, with an additional gradient of a pressure (only living in the saturated zone {ρ = 1}) in the drift. The paper proves existence and some estimates, based on optimal transport techniques.

Introduction

In the past few years modeling crowd behavior has become a very active field of applied mathematics. Beyond their importance in real life applications, these modeling problems serve as basic ideas to understand many other phenomena coming for example from biology (cell migration, tumor growth, pattern formations in animal populations, etc.), particle physics and economics. A first non-exhaustive list of references for these problems is [START_REF] Chalons | Numerical Approximation of a macroscopic model of pedestrian flows[END_REF][START_REF] Colombo | Pedestrian Flows and non-classical shocks[END_REF][START_REF] Coscia | First-Order macroscopic modelling of human crowd dynamics[END_REF][START_REF] Cristiani | Multiscale Modeling of Pedestrian Dynamics[END_REF][START_REF] Dogbe | On the numerical solutions of second order macroscopic models of pedestrian flows[END_REF][START_REF] Helbing | A fluid dynamic model for the movement of pedestrians[END_REF][START_REF] Helbing | Social force model for pedestrian dynamics[END_REF][START_REF] Hughes | A continuum theory for the flow of pedestrian[END_REF][START_REF] Hughes | The flow of human crowds[END_REF][START_REF] Maury | Handling of contacts in crowd motion simulations[END_REF]. A very natural question in all these models is the problem of congestion phenomenon: in many practical situations, very high quantities of individuals could try to occupy the same spot, which could be impossible, or lead to strong negative effects on the motion, because of natural limitations on the crowd density.

These phenomena have been studied by using different models, which could be either "microscopic" (based on ODEs on the motion of a high number of agents) or "macroscopic" (describing the agents via their density and velocity, typically with Eulerian formalism). Let us concentrate on the macroscopic models, where the density ρ plays a crucial role. These very same models can be characterized either by "soft congestion" effects (i.e. the higher the density the slower the motion), or by "hard congestion" (i.e. an abrupt threshold effect: if the density touches a certain maximal value, the motion is strongly affected, while nothing happens for smaller values of the density). See [START_REF] Maury | Handling congestion in crowd motion modeling[END_REF] for comparison between the different classes of models. This last class of models, due to the discontinuity in the congestion effects, presents new mathematical difficulties, which cannot be analyzed with the usual techniques from conservation laws (or, more generally, evolution PDEs) used for soft congestion.

A very powerful tool to attack macroscopic hard-congestion problems is the theory of optimal transportation (see [START_REF] Villani | Topics in Optimal Transportation[END_REF][START_REF] Santambrogio | Optimal Transport for Applied Mathematicians[END_REF]), as we can see in [START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF][START_REF] Maury | Handling congestion in crowd motion modeling[END_REF][START_REF] Roudneff-Chupin | Modélisation macroscopique de mouvements de foule[END_REF][START_REF] Santambrogio | Gradient flows in Wasserstein spaces and applications to crowd movement[END_REF]. In this framework, the density of the agents solves a continuity equation (with velocity field taking into account the congestion effects), and can be seen as a curve in the Wasserstein space.

Our aim in this paper is to endow the macroscopic hard congestion models of [START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF][START_REF] Maury | Handling congestion in crowd motion modeling[END_REF][START_REF] Roudneff-Chupin | Modélisation macroscopique de mouvements de foule[END_REF][START_REF] Santambrogio | Gradient flows in Wasserstein spaces and applications to crowd movement[END_REF]] with diffusion effects. In other words, we will study an evolution equation where particles

• have a spontaneous velocity field u t (x) which depends on time and on their position, and is the velocity they would follow in the absence of the other particles, • must adapt their velocity to the existence of an incompressibility constraint which prevents the density to go beyond a given threshold, • are subject to some diffusion effect. This can be considered as a model for a crowd where a part of the motion of each agent is driven by a Brownian motion. Implementing this new element into the existing models could give a better approximation of reality: as usual when one adds a stochastic component, this can be a (very) rough approximation of unpredictable effects which are not already handled by the model, and this could work well when dealing with large populations.

Anyway, we do not want to discuss here the validity of this hard-congestion model and we are mainly concerned with its mathematical analysis. In particular, we will consider existence and regularity estimates, while we do not treat the uniqueness issue. Uniqueness is considered in a recent work of the first author in collaboration with S. Di Marino, see [START_REF] Marino | Uniqueness issues for evolution equations with density constraints[END_REF], and one can observe that the insertion of diffusion dramatically simplifies the picture as far as uniqueness is concerned.

We also underline that one of the goals of the current paper (and of [START_REF] Marino | Uniqueness issues for evolution equations with density constraints[END_REF]) is to better "prepare" these hard congestion crowd motion models for a possible analysis in the framework of Mean Field Games (see [START_REF] Lasry | Jeux à champ moyen I. Le cas stationnaire[END_REF][START_REF] Lasry | Jeux à champ moyen II. Horizon fini et contrôle optimal[END_REF][START_REF] Lasry | Mean field games[END_REF], and also [START_REF] Santambrogio | A Modest Proposal for MFG with Density Constraints, proceedings of the conference Mean Field Games and related Topics, Roma 1[END_REF]). These MFG models usually involve a stochastic term, also implying regularizing effects, that are useful in the mathematical analysis of the corresponding PDEs.

1.1. The existing first order models in the light of [START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF][START_REF] Maury | Handling congestion in crowd motion modeling[END_REF]. Some macroscopic models for crowd motion with density constraints and "hard congestion" effects were studied in [START_REF] Maury | Handling congestion in crowd motion modeling[END_REF] and [START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF]. We briefly present them as follows:

• The density of the population in a bounded (convex) domain Ω ⊂ R d is described by a probability measure ρ ∈ P(Ω). The initial density ρ 0 ∈ P(Ω) evolves in time, and ρ t denotes its value at each time t ∈ [0, T ]. • The spontaneous velocity field of the population is a given time-dependent field, denoted by u t . It represents the velocity that each individual would like to follow in the absence of the others. Ignoring the density constraint, this would give rise to the continuity equation ∂ t ρ t + ∇ • (ρ t u t ) = 0. We observe that in the original work [START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF] the vector field u t (x) was taken of the form -∇D(x) (independent of time and of gradient form) but we try here to be more general (see [START_REF] Roudneff-Chupin | Modélisation macroscopique de mouvements de foule[END_REF] where the non-gradient case is studied under some stronger regularity assumptions). • The set of admissible densities will be denoted by K := {ρ ∈ P(Ω) : ρ ≤ 1}. In order to guarantee that K is neither empty nor trivial, we suppose |Ω| > 1.

• The set of admissible velocity fields with respect to the density ρ is characterized by the sign of the divergence of the velocity field on the saturated zone. We need to suppose also that all admissible velocity fields are such that no mass exists from the domain. So formally we set adm(ρ) := v : Ω → R d : ∇ • v ≥ 0 on {ρ = 1} and v • n ≤ 0 on ∂Ω .

• We consider the projection operator P in L 2 (L d ):

P adm(ρ) [u] ∈ argmin v∈adm(ρ) Ω |u -v| 2 dx.
Note that we could have used the Hilbert space L 2 (ρ) instead of L 2 (L d ): this would be more natural in this kind of evolution equations, as L 2 (ρ) is interpreted in a standard way as the tangent space to the Wasserstein space W 2 (Ω). Yet, these two projections turn out to be the same in this case, as the only relevant zone is {ρ = 1}. This is just formal, and would require more rigorous definitions (in particular of the divergence constraint in adm(ρ), see below). Anyway, to clarify, we choose to use the L 2 (L d )-projection: in this way the vector fields are considered as defined Lebesgue-a.e. on the whole Ω (and not only on {ρ > 0}) and the dependence of the projected vector field on ρ only passes through the set adm(ρ). • Finally we solve the following modified continuity equation for ρ

(1.1) ∂ t ρ t + ∇ • ρ t P adm(ρt) [u t ] = 0,
where the main point is that ρ is advected by a vector field, compatible with the constraints, which is the closest to the spontaneous one.

The problem in solving Equation (1.1) is that the projected field has very low regularity: it is a priori only L 2 in x, and it does not depend smoothly on ρ either (since a density 1 and a density 1 -ε give very different projection operators). By the way, its divergence is not well-defined either.

To handle this issue we need to redefine the set of admissible velocities by duality. Taking a test function p ∈ H 1 (Ω), p ≥ 0 a.e., we obtain by the integration-by-parts equality

Ω v • ∇p dx = - Ω (∇ • v)p dx + ∂Ω pv • n dH d-1 (x).
For vector fields v which do not let mass go through the boundary ∂Ω we have (in an a.e. sense) v • n = 0. This leads to the following definition

adm(ρ) = v ∈ L 2 (Ω; R d ) : Ω v • ∇p dx ≤ 0, ∀p ∈ H 1 (Ω), p ≥ 0, p(1 -ρ) = 0 a.e. ,
(indeed, for smooth vector field with vanishing normal component on the boundary, this is equivalent to imposing ∇ • v ≥ 0 on the set {ρ = 1}). Now, if we set

press(ρ) := p ∈ H 1 (Ω) : p ≥ 0, p(1 -ρ) = 0 a.e. ,
we observe that, by definition, adm(ρ) and ∇press(ρ) are two convex cones which are dual to each other in L 2 (Ω; R d ). Hence we always have a unique orthogonal decomposition

(1.2) u = v + ∇p, v ∈ adm(ρ); p ∈ press(ρ), Ω v • ∇p dx = 0.
In this decomposition (as it is the case every time we decompose on two dual convex cones), v = P adm(ρ) [u]. These will be our mathematical definitions for adm(ρ) and for the projection onto this cone. Via this approach (introducing the new variable p and using its characterization from the previous line), for a given desired velocity field u : [0, T ] × Ω → R d , the continuity equation (1.1) can be rewritten as a system for the pair of variables (ρ, p) which is

(1.3)    ∂ t ρ t + ∇ • (ρ t (u t -∇p t )) = 0, in [0, T ] × Ω, p ≥ 0, ρ ≤ 1, p(1 -ρ) = 0, in [0, T ] × Ω, ρ t (u t -∇p t ) • n = 0, on [0, T ] × ∂Ω.
This system is endowed with the initial condition ρ(0, x) = ρ 0 (x) (ρ 0 ∈ K). As far as the spatial boundary ∂Ω is concerned, we put no-flux boundary conditions to preserve the mass in Ω. Note that in the above system we withdrew the condition (u t -∇p t ) • ∇p t = 0, as it is a consequence of the system (1.3) itself. Informally, this can be seen in the following way: for an arbitrary p 0 ∈ press(ρ t 0 ), we have that t → Ω p 0 ρ t is maximal at t = t 0 (where it is equal to Ω p 0 ). Differentiating this quantity w.r.t. t at t = t 0 , using the equation (1.3), we get the desired orthogonality condition at t = t 0 . For a rigorous proof of this fact (which holds for a.e. t 0 ), we refer to Proposition 4.7 in [START_REF] Marino | Measure sweeping processes[END_REF]. 1.2. A diffusive counterpart. The goal of our work is to study a second order model of crowd movements with hard congestion effects where beside the transport factor a non-degenerate diffusion is present as well. The diffusion is the consequence of a randomness (a Brownian motion) in the movement of the crowd.

With the ingredients that we introduced so far, we will modify the Fokker-Planck equation ∂ t ρ t -∆ρ t + ∇ • (ρ t u t ) = 0 in order to take into account the density constraint ρ t ≤ 1. Assuming enough regularity for the velocity field u, we observe that the Fokker Planck equation is derived from a motion given by the SODE dX t = u t (X t ) dt + √ 2 dB t (where B t is the standard d-dimensional Brownian motion), but is macroscopically represented by the advection of the density ρ t by the vector field -∇ρ t /ρ t + u t . Projecting onto the set of admissible velocities raises a natural question: should we project only u t , and then apply the diffusion, or project the whole vector field, including -∇ρ t /ρ t ? But this is not a real issue, since, at least formally, ∇ρ t /ρ t = 0 on the saturated set {ρ t = 1} and P adm(ρt) [-∇ρ t /ρ t + u t ] = P adm(ρt) [-∇ρ t /ρ t ] + P adm(ρt) [u t ] = 0 + P adm(ρt) [u t ]. Rigorously, this corresponds to the fact that the Heat Kernel preserves the constraint ρ ≤ 1. As a consequence, we consider the modified Fokker-Planck type equation

(1.4) ∂ t ρ t -∆ρ t + ∇ • ρ t P adm(ρt) [u t ] = 0,
which can also be written equivalently for the variables (ρ, p) as

(1.5) ∂ t ρ t -∆ρ t + ∇ • (ρ t (u t -∇p t )) = 0, in [0, T ] × Ω, p ≥ 0, ρ ≤ 1, p(1 -ρ) = 0, in [0, T ] × Ω.
As usual, these equations are complemented by no-flux boundary conditions and by an initial datum ρ(0, x) = ρ 0 (x).

Roughly speaking, we can consider that this equation describes the law of a motion where each agent solves the stochastic differential equation

dX t = (u t (X t ) -∇p t (X t )) dt + √ 2 dB t .
This last statement is just formal and there are several issues defining an SODE like this: indeed, the pressure variable is also an unknown, and globally depends on the law ρ t of X t . Hence, if we wanted to see this evolution as a superposition of individual motions, each agent should somehow predict the evolution of the pressure in order to solve his own equation. This reminds of some notions from the stochastic control formulation of Mean-Field Games, as introduced by J.-M. Lasry and P.-L.

Lions, even if here there are no strategic issues for the players. For MFG with density constraints, we refer to [START_REF] Cardaliaguet | First order Mean Field Games with density constraints: Pressure equals Price[END_REF][START_REF] Mészáros | A variational approach to second order Mean Field Games with density constraints: the stationary case[END_REF][START_REF] Santambrogio | A Modest Proposal for MFG with Density Constraints, proceedings of the conference Mean Field Games and related Topics, Roma 1[END_REF]. However, in this paper we will not consider any microscopic or individual problem, but only study the parabolic PDE (1.5).

1.3. Structure of the paper and main results. The main goal of the paper is to provide an existence result, with some extra estimates, for the Fokker-Planck equation (1.5) via time discretization, using the so-called splitting method (the two main ingredients of the equation, i.e. the advection with diffusion on one hand, and the density constraint on the other hand, are treated one after the other). In Section 2 we will collect some preliminary results, including what we need from optimal transport and from the previous works about density-constrained crowd motion, in particular on the projection operator onto the set K. In Section 3 we will provide the existence result we aim at, by a splitting scheme and some entropy bounds; the solution will be a curve of measures in AC 2 ([0, T ]; W 2 (Ω)) (absolutely continuous curves with square-integrable speed). In Section 4 we will make use of BV estimates to justify that the solution we just built is also Lip([0, T ]; W 1 (Ω)) and satisfies a global BV bound ρ t BV ≤ C (provided ρ 0 ∈ BV ): this requires to combine BV estimates on the Fokker-Planck equation (which are available depending on the regularity of the vector field u) with BV estimates on the projection operator on K (which have been recently proven in [START_REF] De Philippis | BV estimates in optimal transportation and applications[END_REF]). Section 5 presents a short review of alternative approaches, all discretized in time, but based either on gradient-flow techniques (the JKO scheme, see [START_REF] Jordan | The variational formulation of the Fokker-Plack equation[END_REF]) or on different splitting methods. Finally, in the Appendix A we detail the BV estimates on the Fokker-Planck equation (without any density constraint) that we could find; this seems to be a delicate matter, interesting in itself, and we are not aware of the sharp assumptions on the vector field u to guarantee the BV estimate that we need.

Preliminaries

2.1. Basic definitions and general facts on optimal transport. Here we collect some tools from the theory of optimal transportation, Wasserstein spaces, its dynamical formulation, etc. which will be used later on. We set our problem either in a compact convex domain Ω ⊂ R d with smooth boundary or in the d-dimensional flat torus Ω := T d (even if we will not adapt all our notations to the torus case). We refer to [START_REF] Villani | Topics in Optimal Transportation[END_REF][START_REF] Santambrogio | Optimal Transport for Applied Mathematicians[END_REF] for more details. Given two probability measures µ, ν ∈ P(Ω) and for p ≥ 1 we define the usual Wasserstein metric by means of the Monge-Kantorovich optimal transportation problem

W p (µ, ν) := inf Ω×Ω |x -y| p dγ(x, y) : γ ∈ Π(µ, ν) 1 p
, where Π(µ, ν) := {γ ∈ P(Ω × Ω) : (π x ) # γ = µ, (π y ) # γ = ν} and π x and π y denote the canonical projections from Ω × Ω onto Ω. This quantity happens to be a distance on P(Ω) which metrizes the weak- * convergence of probability measures; we denote by W p (Ω) := (P(Ω), W p ), i.e. the space of probabilities on Ω endowed with this distance.

Moreover, in the quadratic case p = 2 and under the assumption µ ≪ L d (the d-dimensional Lebesgue measure on Ω) in the late 80's Y. Brenier showed (see [START_REF] Brenier | Décomposition polaire et réarrangement monotone des champs de vecteurs[END_REF][START_REF] Brenier | Polar factorization and monotone rearrangement of vector valued functions[END_REF]) that actually the optimal γ in the above problem (the existence of which is obtained simply by the direct method of calculus of variations) is induced by a map, which is the gradient of a convex function, i.e. there exists S : Ω → Ω and ψ : Ω → R convex such that S = ∇ψ and γ := (id, S) # µ. The function ψ is obtained as ψ(x) = 1 2 |x| 2 -ϕ(x), where ϕ is the so-called Kantorovich potential for the transport from µ to ν, and is characterized as the solution of a dual problem that we will not develop here. In this way, the optimal transport map S can also be written as S(x) = x -∇ϕ(x). Later in the 90's R. McCann (see [START_REF] Mccann | A convexity principle for interacting gases[END_REF]) introduced a notion of interpolation between probability measures: the curve µ t := ((T -t)x + ty) # γ, for t ∈ [0, T ] (T > 0 is given), gives a constant speed geodesic in the Wasserstein space connecting µ 0 := µ and µ T := ν.

Based on this notion of interpolation in 2000 J.-D. Benamou and Y. Brenier used some ideas from fluid mechanics to give a dynamical formulation to the Monge-Kantorovich problem (see [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF]). They showed that

1 pT p-1 W p p (µ, ν) = inf {B p (E, µ) : ∂ t µ + ∇ • E = 0, µ 0 = µ, µ T = ν} .
Here B p is a functional defined on pairs (E, µ), where E is a d-dimensional vector measure on [0, T ] × Ω and µ = (µ t ) t is a Borel-measurable family of probability measures on Ω. This functional is defined to be finite only if E = E t ⊗ dt (i.e. it is induced by a measurable family of vector measures on Ω:

[0,T ]×Ω ξ(t, x)• dE(t, x) = T 0 dt Ω ξ(t, x)• dE t (x) for all test functions ξ ∈ C 0 ([0, T ]×Ω; R d ))
and in this case it is defined through

B p (E, µ) :=    T 0 Ω 1 p |v t | p dµ t (x) dt, if E t = v t • µ t +∞, otherwise.
It is well-known that B p is jointly convex and l.s.c. w.r.t the weak- * convergence of measures (see Section 5.3.1 in [START_REF] Santambrogio | Optimal Transport for Applied Mathematicians[END_REF]) and that, if

∂ t µ+∇•E = 0, then B p (E, µ) < +∞ implies that t → µ t is a curve in AC p ([0, T ]; W p (Ω)) 1 .
In particular it is a continuous curve and the initial and final conditions on µ 0 and µ T are well-defined.

Coming back to curves in Wasserstein spaces, it is well known (see [START_REF] Ambrosio | Gradient flows in metric spaces and in the space of probability measure[END_REF] or Section 5.3 in [START_REF] Santambrogio | Optimal Transport for Applied Mathematicians[END_REF]) that for any distributional solution µ t (being a continuous curve in W p (Ω)) of the continuity equation

∂ t µ + ∇ • E = 0 with E t = v t • µ t , we have the relations |µ ′ | Wp (t) ≤ v t L p µ t and W p (µ t , µ s ) ≤ t s |µ ′ | Wp (τ ) dτ,
where we denoted by |µ ′ | Wp (t) the metric derivative w.r.t. W p of the curve µ t (see for instance [START_REF] Ambrosio | Topics on analysis in metric spaces[END_REF] for general notions about curves in metric spaces and their metric derivative). For curves µ t that are geodesics in W p (Ω) we have the equality

W p (µ 0 , µ 1 ) = 1 0 |µ ′ | Wp (t) dt = 1 0 v t L p µ t dt.
The last equality is in fact the Benamou-Brenier formula with the optimal velocity field v t being the density of the optimal E t w.r.t. the optimal µ t . This optimal velocity field v t can be computed as v t := (S -id) • (S t ) -1 , where S t := (1 -t)id + tS is the transport in McCann's interpolation (we assume here that the initial measure µ 0 is absolutely continuous, so that we can use transport maps instead of plans). This expression can be obtained if we consider that in this interpolation particles move with constant speed S(x) -x, but x represents here a Lagrangian coordinate, and not an Eulerian one: if we want to know the velocity at time t at a given point, we have to find out first the original position of the particle passing through that point at that time. 1 Here AC p ([0, T ]; Wp(Ω)) denotes the class of absolutely continuous curves in Wp(Ω) with metric derivative in L p . See the connection with the functional Bp.

In the sequel we will also need the notion of entropy of a probability density, and for any probability measure ̺ ∈ P(Ω) we define it as

E(̺) :=    Ω ̺(x) log ̺(x) dx, if ̺ ≪ L d , +∞, otherwise.
We recall that this functional is l.s.c. and geodesically convex in W 2 (Ω). As we will be mainly working with absolutely continuous probability measures (w.r.t. Lebesgue), we often identify measures with their densities.

2.2. Projection problems in Wasserstein spaces. Our analysis strongly relies on the projection operator P K in the sense of W 2 . Here K := {ρ ∈ P(Ω) : ρ ≤ 1} and

P K [µ] := argmin ρ∈K 1 2 W 2 2 (µ, ρ).
We recall (see [START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF][START_REF] Santambrogio | Gradient flows in Wasserstein spaces and applications to crowd movement[END_REF] and [START_REF] De Philippis | BV estimates in optimal transportation and applications[END_REF]) the main properties of the projection P K operator.

• As far as Ω is compact, for any probability measure µ, the minimizer in min ρ∈K 1 2 W 2 2 (µ, ρ) exists and is unique, and the operator P K is continuous (it is even C 0,1/2 for the W 2 distance).

• The projection P K [µ] saturates the constraint ρ ≤ 1 in the sense that for any µ ∈ P(Ω) there exists a measurable set B ⊆ Ω such that

P K [µ] = ½ B + µ ac ½ B c
, where µ ac is the absolutely continuous part of µ. • The projection is characterized in terms of a pressure field, in the sense that ρ = P K [µ] if and only if there exists a Lipschitz function p ≥ 0, with p(1 -ρ) = 0, and such that the optimal transport map S from ρ to µ is given by S := id -∇ϕ = id + ∇p. • There is (as proven in [START_REF] De Philippis | BV estimates in optimal transportation and applications[END_REF]) a quantified BV estimate: if µ ∈ BV (in the sense that it is absolutely continuous and that its density belongs to BV (Ω)), then P K [µ] is also BV and

T V (P K [µ], Ω) ≤ T V (µ, Ω).
This last BV estimate will be crucial in Section 4, and it is important to have it in this very form (other estimates of the form T V (P K [µ], Ω) ≤ aT V (µ, Ω) + b would not be as useful as this one, as they cannot be easily iterated).

Existence via a splitting-up type algorithm (Main Scheme)

Similarly to the approach in [START_REF] Maury | Handling congestion in crowd motion modeling[END_REF] (see the algorithm (13) and Theorem 3.5) for a general, nongradient, vector field, we will build a theoretical algorithm, after time-discretization, to produce a solution of (1.5). Let us remark that splitting-type methods have been widely used in other contexts as well, see for instance the paper [START_REF] Ph | A Trotter product formula for gradient flows in metric spaces[END_REF] which deals with splitting methods for Fokker-Planck equations and for more general gradient flows in metric and Wasserstein spaces, or [START_REF] Laborde | On some non linear evolution systems which are perturbations of Wasserstein gradient flows[END_REF] where a splitting-like approach is used to attack PDEs which are not gradient flows but "perturbations" of gradient flows.

In this section the spontaneous velocity field is a general vector field u : [0, T ] × Ω → R d (not necessarily a gradient), which depends also on time. The only assumption we require on u is the following:

(U) u ∈ L ∞ ([0, T ] × Ω; R d ).
We will work on a time interval [0, T ] and in a bounded convex domain Ω ⊂ R d (the case of the flat torus is even simpler and we will not discuss it in details). We consider ρ 0 ∈ P ac (Ω) to be given, which represents the initial density of the population, and we suppose ρ 0 ∈ K.

3.1. Splitting using the Fokker-Planck equation. Let us consider the following scheme.

Main scheme: Let τ > 0 be a small time step with N := ⌊T /τ ⌋. Let us set ρ τ 0 := ρ 0 and for every k ∈ {1, . . . , N } we define ρ τ k+1 from ρ τ k in the following way. First we solve (3.1)

∂ t ̺ t -∆̺ t + ∇ • (̺ t u t+kτ ) = 0, t ∈]0, τ ], ̺ 0 = ρ τ k , equipped with the no-flux boundary condition (̺ t (∇̺ t -u t ) • n = 0 a.e. on ∂Ω) and set ρ τ k+1 = P K [ρ τ k+1 ]
, where ρτ k+1 = ̺ τ . See Figure 1 on the right.

• ρ τ k ∂ t ̺ t - ∆ ̺ t + ∇ • ( ̺ t u t + k τ ) = 0 • ρτ k+1 = ̺ τ id + τ ∇ p τ k + 1 • ρ τ k+1 Figure 1. One time step
Let us remark first that by classical results on parabolic equations (see for instance [START_REF] Ladyzenskaja | Linear and quasilinear equations of parabolic type, Translated from the Russian[END_REF]), since u satisfies the assumption (U), Problem 3.1 admits a unique distributional solution.

The above algorithm means the following: first follow the Fokker-Planck equation, ignoring the density constraint, for a time τ , then project. In order to state and prove the convergence of the scheme, we need to define some suitable interpolations of the discrete sequence of densities that we have just introduced.

First interpolation. We define the following curves of densities, velocities and momenta constructed with the help of the ρ τ k 's. First set

ρ τ t := ̺ 2(t-kτ ) , if t ∈ [kτ, (k + 1/2)τ [ , id + 2((k + 1)τ -t)∇p τ k+1 # ρ τ k+1 , if t ∈ [(k + 1/2)τ, (k + 1)τ [ ,
where ̺ t is the solution of the Fokker-Planck equation (3.1) with initial datum ρ τ k and ∇p τ k+1 arises from the projection of ρτ k+1 , more precisely (id + τ ∇p τ k+1 ) is the optimal transport from ρ τ k+1 to ρτ k+1 . What are we doing? We are fitting into a time interval of length τ the two steps of our algorithm. First we follow the FP equation (3.1) at double speed, then we interpolate between the measure we reached and its projection following the geodesic between them. This geodesic is easily described as an image measure of ρ τ k+1 through McCann's interpolation. By the construction it is clear that ρ τ t is a continuous curve in P(Ω) for t ∈ [0, T ]. We now define a family of time-dependent vector fields though 2) in the first construction we have taken into account the projection operator explicitly, while in the second one we see it just in an indirect manner (via the 'jumps' occurring at every time of the form t = kτ ). The third interpolation is piece-wise constant, and at every time it satisfies the density constraint;

v τ t := -2 ∇̺ 2(t-kτ ) ̺ 2(t-kτ ) + 2u t , if t ∈ [kτ, (k + 1/2)τ [ , -2∇p τ k+1 • (id + 2((k + 1)τ -t)∇p τ k+1 ) -1 , if t ∈ [(k + 1/2)τ, (k + 1)
3) in the first interpolation the pair (ρ τ , E τ ) solves the continuity equation, while in the other two they do not. This is not astonishing, as the continuity equation characterizes continuous curves in W 2 (Ω).

In order to prove the convergence of the scheme above, we will obtain uniform AC 2 ([0, T ]; W 2 (Ω)) bounds for the curves ρ τ . A key observation here is that the metric derivative (w.r.t. W 2 ) of the solution of the Fokker-Planck equation is comparable with the time differential of the entropy functional along the same solution (see Lemma 3.2). Now we state the main theorem of this section.

Theorem 3.1. Let ρ 0 ∈ K and u be a given desired velocity field satisfying (U). Let us consider the interpolations introduced above. Then there exists a continuous curve [0, T ] ∋ t → ρ t ∈ W 2 (Ω) and some vector measures E, Ẽ, Ê ∈ M([0, T ] × Ω) such that the curves ρ τ , ρτ , ρτ converge uniformly in W 2 (Ω) to ρ and

E τ * ⇀ E, Ẽτ * ⇀ Ẽ, Êτ * ⇀ Ê, in M([0, T ] × Ω) d , as τ → 0.
Moreover E = Ẽ -Ê and for a.e. t ∈ [0, T ] there exist time-dependent measurable vector fields v t , ṽt , vt such that

(1) E = ρv, Ẽ = ρṽ, Ê = ρv, (2)
T 0 v t 2 L 2 ρ t + ṽt 2 L 2 ρ t + vt 2 L 2 ρ t dt < +∞, (3)
v t = ṽt -vt ρ t -a.e., Ẽt = ρ t u t -∇ρ t and vt = ∇p t , ρ t -a.e., where p ∈ L 2 ([0, T ]; H 1 (Ω)), p ≥ 0 and p(1 -ρ) = 0 a.e. in [0, T ] × Ω. As a consequence, the pair (ρ, p) is a weak solution of the problem

(3.2)        ∂ t ρ t -∆ρ t + ∇ • (ρ t (u t -∇p t )) = 0, in [0, T ] × Ω, p t ≥ 0, ρ t ≤ 1, p t (1 -ρ t ) = 0, in [0, T ] × Ω, ρ t (∇ρ t -u t + ∇p t ) • n = 0, on [0, T ] × ∂Ω, ρ(0, •) = ρ 0 .
To prove this theorem we will use the following tools.

(3.3) 1 2 b a Ω - ∇̺ t ̺ t + u t 2 ̺ t dx dt ≤ E(̺ a ) -E(̺ b ) + 1 2 b a Ω |u t | 2 ̺ t dx dt
In particular this implies

(3.4) 1 2 b a |̺ ′ t | 2 W 2 dt ≤ E(̺ a ) -E(̺ b ) + 1 2 b a Ω |u t | 2 ̺ t dx dt,
where

|̺ ′ t | W 2 denotes the metric derivative of the curve t → ̺ t ∈ W 2 (Ω). Proof.
To prove this inequality, we will first make computations in the case where both u and ̺ are smooth, and ̺ is boundeded from below by a positive constant. In this case we can write

d dt E(̺ t ) = Ω (log ̺ t + 1)∂ t ̺ t dx = Ω log ̺ t (∆̺ t -∇ • (̺ t u t )) dx = Ω - |∇̺ t | 2 ̺ t + u t • ∇̺ t dx,
where we used the conservation of mass (i.e. Ω ∂ t ̺ t dx = 0) and the boundary conditions in the integration by parts. We now compare this with

1 2 Ω - ∇̺ t ̺ t + u t 2 ̺ t dx - 1 2 Ω |u t | 2 ̺ t dx = Ω 1 2 |∇̺ t | 2 ̺ t -∇̺ t • u t dx ≤ Ω |∇̺ t | 2 ̺ t -∇̺ t • u t dx = - d dt E(̺ t ).
This provides the first part of the statement, i.e. (3.3). If we combine this with the fact that the metric derivative of the curve t → ̺ t is always less or equal than the L 2 ̺t norm of the velocity field in the continuity equation, we also get 1 2

|̺ ′ t | 2 W 2 - 1 2 Ω |u t | 2 ̺ t ≤ - d dt E(̺ t ),
and hence (3.4).

In order to prove the same estimates without artificial smoothness and lower bound assumptions, we can act by approximation. We approximate the density ̺ a by smooth and strictly positive densities ̺ k a (by convolution, so that we guarantee in particular E(̺ k a ) → E(̺ a )), and the vector field u with smooth vector fields u k (strongly in L

E(̺ b ) -E(̺ a ) ≤ 1 2 b a Ω |u t | 2 ̺ t dx dt,
hence in particular for u satisfying (U) we have

E(̺ b ) -E(̺ a ) ≤ 1 2 u 2 L ∞ (b -a). As a consequence, if ̺ a ≤ 1, then we have E(̺ b ) ≤ 1 2 u 2 L ∞ (b -a).
The same estimate can be applied to the curve ρτ , with a = kτ and b ∈]kτ, (k + 1)τ [, thus obtaining E(ρ τ t ) ≤ Cτ for every t. Lemma 3.4. For any ρ ∈ P(Ω) we have

E (P K [ρ]) ≤ E(ρ).
Proof. We can assume ρ ≪ L d , otherwise the claim is straightforward. As we pointed out in Section 2.2, we know that there exists a measurable set B ⊆ Ω such that

P K [ρ] = ½ B + ρ½ B c .
Hence it is enough to prove that

B ρ log ρ dx ≥ 0 = B P K [ρ] log P K [ρ] dx,
as the entropies on B c coincide. As the mass of ρ and P K [ρ] is the same on the whole Ω, and they coincide on B c , we have

B ρ(x) dx = B P K [ρ] dx = |B|.
Then, by Jensen's inequality we have 1

|B| B ρ log ρ dx ≥ 1 |B| B ρ dx log 1 |B| B ρ dx = 0.
The entropy decay follows.

To analyze the pressure field we will need the following result.

Lemma 3.5. Let {p τ } τ >0 be a bounded sequence in L 2 ([0, T ]; H 1 (Ω)) and {ρ τ } τ >0 a sequence of piecewise constant curves valued in W 2 (Ω), which satisfy W 2 (ρ τ (a), ρ τ (b)) ≤ C √ b -a + τ for all a < b ∈ [0, T ] for a fixed constant C. Suppose that p τ ≥ 0, p τ (1 -ρ τ ) = 0, ρ τ ≤ 1,
and that

p τ ⇀ p weakly in L 2 ([0, T ]; H 1 (Ω)) and ρ τ → ρ uniformly in W 2 (Ω). Then p(1 -ρ) = 0 a.e. in [0, T ] × Ω.
The proof of this result is the same as in Step 3 of Section 3.2 of [START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF] (see also [START_REF] Roudneff-Chupin | Modélisation macroscopique de mouvements de foule[END_REF] and Lemma 4.6 in [START_REF] Marino | Measure sweeping processes[END_REF]). We omit it in order not to overburden the paper.

The reader can note the strong connection with the classical Aubin-Lions lemma [START_REF] Aubin | Un théorème de compacité[END_REF], applied to the compact injection of L 2 into H -1 . Indeed, from the weak convergence of p τ in L 2 ([0, T ]; H 1 (Ω)), we just need to provide strong convergence of ρ τ in L 2 ([0, T ]; H -1 (Ω)). If instead of the quasi-Hölder assumption of the above lemma we suppose a uniform bound of {ρ τ } τ in AC 2 ([0, T ]; W 2 (Ω)) (which is not so different), then the statement can be really deduced from the Aubin-Lions lemma. Indeed, the sequence {ρ τ } is bounded in L ∞ ([0, T ]; L 2 (Ω)) and its time-derivative would be bounded in L 2 ([0, T ]; H -1 (Ω)). This strongly depends on the fact that the H -1 distance can be controlled by the W 2 distance as soon as the measures have uniformly bounded densities (see [START_REF] Loeper | Uniqueness of the solution to the Vlasov-Poisson system with bounded density[END_REF][START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF]), a tool which also crucial in the proofs in [START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF][START_REF] Roudneff-Chupin | Modélisation macroscopique de mouvements de foule[END_REF][START_REF] Marino | Measure sweeping processes[END_REF]. Then, the Aubin-Lions lemma guarantees compactness in C 0 ([0, T ]; H -1 (Ω)), which is more than what we need. Lemma 3.6. Let us consider the previously defined interpolations. Then we have the following facts.

(i) For every τ > 0 and k we have

max W 2 2 (ρ τ k , ρτ k+1 ), W 2 2 (ρ τ k , ρ τ k+1 ) ≤ τ C E(ρ τ k ) -E(ρ τ k+1 ) + Cτ 2 , where C > 0 only depends on u L ∞ . (ii) There exists a constant C, only depending on ρ 0 and u L ∞ , such that B 2 (E τ , ρ τ ) ≤ C, B 2 ( Ẽτ , ρτ ) ≤ C and B 2 ( Êτ , ρτ ) ≤ C. (iii) For the curve [0, T ] ∋ t → ρ τ t we have that T 0 |(ρ τ t ) ′ | 2 W 2 dt ≤ C,
= P K [ρ τ k+1 ] we have that (3.5) W 2 (ρ τ k , ρ τ k+1 ) ≤ W 2 (ρ τ k , ρτ k+1 ) + W 2 (ρ τ k+1 , ρ τ k+1 ) ≤ 2W 2 (ρ τ k , ρτ k+1
). We use (as before) the notation ̺ t , t ∈ [0, τ ] for the solution of the Fokker-Planck equation (3.1) with initial datum ρ τ k , in particular we have ̺ τ = ρτ k+1 . Using Lemma 3.2 and since ̺ 0 = ρ τ k and ̺ τ = ρτ k+1 we have by (3.4) and using

W 2 (ρ τ k , ρτ k+1 ) ≤ τ 0 |̺ ′ t | W 2 dt W 2 2 (ρ τ k , ρτ k+1 ) ≤ τ 1 2 τ 0 |̺ ′ t | 2 W 2 dt 1 2 2 ≤ 2τ (E(̺ 0 ) -E(̺ τ )) + τ τ 0 Ω |u kτ +t | 2 ̺ t dx dt ≤ 2τ E(ρ τ k ) -E(ρ τ k+1 ) + Cτ 2 ≤ 2τ E(ρ τ k ) -E(ρ τ k+1 ) + Cτ 2
, where C > 0 is a constant depending just on u L ∞ . We also used the fact that E(ρ τ k+1 ) ≤ E(ρ τ k+1 ), a consequence of Lemma 3.4. Now by the means of (3.5) we obtain

W 2 2 (ρ τ k , ρ τ k+1 ) ≤ τ C E(ρ τ k ) -E(ρ τ k+1 ) + Cτ 2 .
(ii) We use Lemma 3.2 on the intervals of type [kτ, (k + 1/2)τ [ and the fact that on each interval of type [(k + 1/2)τ, (k + 1)τ [ the curve ρ τ t is a constant speed geodesic. In particular, on these intervals we have

|(ρ τ ) ′ | W 2 = v τ t L 2 ρ τ t = 2τ ∇p τ k+1 L 2 ρ τ k+1 = 2W 2 (ρ τ k+1 , ρτ k+1 ).
On the other hand we also have

τ 2 ∇p τ k+1 2 L 2 ρ τ k+1 = W 2 2 (ρ τ k+1 , ρτ k+1 ) ≤ W 2 2 (ρ τ k , ρτ k+1 ) ≤ τ C E(ρ τ k ) -E(ρ τ k+1 ) + Cτ 2 .
Hence we obtain

(k+1)τ kτ v τ t 2 L 2 (ρ τ t ) dt = (k+1/2)τ kτ Ω 4 - ∇̺ 2(t-kτ ) ̺ 2(t-kτ ) + u 2t-kτ 2 ̺ 2(t-kτ ) (x) dx dt + 4 (k+1)τ (k+1/2)τ Ω |∇p τ k+1 | 2 ρ τ k+1 dx dt ≤ C E(ρ τ k ) -E(ρ τ k+1 ) + Cτ + 2τ ∇p τ k+1 2 L 2 ρ τ k+1 ≤ C E(ρ τ k ) -E(ρ τ k+1 ) + Cτ.
Hence by adding up we obtain

B 2 (E τ , ρ τ ) ≤ k C E(ρ τ k ) -E(ρ τ k+1 ) + Cτ = C E(ρ τ 0 ) -E(ρ τ N +1 ) + CT ≤ C.
The estimate on B 2 ( Ẽτ , ρτ ) and B 2 ( Êτ , ρτ ) are completely analogous and descend from the previous computations.

(iii) The estimate on B 2 (E τ , ρ τ ) implies a bound on

T 0 |(ρ τ t ) ′ | 2
W 2 dt because v τ is a velocity field for ρ τ (i.e., the pair (E τ , ρ τ ) solves the continuity equation).

(iv) In order to estimate the total mass of E we write

|E τ |([0, T ] × Ω) = T 0 Ω |v τ t |ρ τ t dx dt ≤ T 0 Ω |v τ t | 2 ρ τ t dx 1 2 Ω ρ τ t dx 1 2 dt ≤ √ T T 0 Ω |v τ t | 2 ρ τ t dx dt 1 2 ≤ C.
The bounds on Ẽτ and Êτ rely on the same argument.

Proof of Theorem 3.1. We use the tools from Lemma 3.6.

Step 1. By the bounds on the metric derivative of the curves ρ τ t we get compactness, i.e. there exists a curve [0, T ] ∋ t → ρ t ∈ P(Ω) such that ρ τ (up to subsequences) converges uniformly in [0, T ] w.r.t. W 2 , in particular weakly- * in P(Ω) for all t ∈ [0, T ]. It is easy to see that ρτ and ρτ are converging to the same curve. Indeed we have ρτ t = ρ τ s(t) and ρτ t = ρ τ ŝ(t) for |s(t) -t| ≤ τ and

|ŝ(t) -t| ≤ τ , which implies W 2 (ρ τ t , ρτ t ), W 2 (ρ τ t , ρτ t ) ≤ Cτ 1 2
. This provides the convergence to the same limit.

Step 2. By the boundedness of E τ , Ẽτ and Êτ in M([0, T ] × Ω) d we have the existence of

E, Ẽ, Ê ∈ M([0, T ] × Ω) d such that (up to a subsequence) E τ * ⇀ E, Ẽτ * ⇀ Ẽ, Êτ * ⇀ Ê as τ → 0.
Now we show that E = Ẽ -Ê. Indeed, let us show that for any test function

f ∈ Lip([0, T ] × Ω) d we have T 0 Ω f t • E τ t -( Ẽτ t + Êτ t ) ( dx, dt) → 0,
as τ → 0. First for each k ∈ {0, . . . , N } we have that

(k+1/2)τ kτ Ω f t • E τ t ( dx, dt) = (k+1)τ kτ Ω f (t+kτ )/2 • (-∇̺ t-kτ + u t ̺ t-kτ )( dx, dt) = (k+1)τ kτ Ω f t • Ẽτ t ( dx, dt) + (k+1)τ kτ Ω f (t+kτ )/2 -f t • Ẽτ t ( dx, dt)
and

(k+1)τ (k+1/2)τ Ω f t • E τ t ( dx, dt) = (k+1)τ kτ Ω -f (t+(k+1)τ )/2 • (id + ((k + 1)τ -t)∇p τ k+1 ) • ∇p τ k+1 ρ τ k+1 ( dx, dt) = - (k+1)τ kτ Ω f t • Êτ t ( dx, dt) + (k+1)τ kτ Ω f t -f (t+(k+1)τ )/2 • (id + ((k + 1)τ -t)) • vτ t ρτ t ( dx, dt)
This implies that

T 0 Ω f t • (E τ t -Ẽτ t + Êτ t )( dx, dt) ≤ k (k+1)τ kτ Lip(f )τ Ω | Ẽτ t |( dx, dt) + k (k+1)τ kτ Lip(f )τ Ω (1 + |v τ t |)| Êτ t |( dx, dt) ≤ τ CLip(f ) | Ẽτ |([0, T ] × Ω) + | Êτ |([0, T ] × Ω) + B 2 ( Ê, ρ) ≤ τ CLip(f ),
for a uniform constant C > 0. Letting τ → 0 we prove the claim.

Step 3. The bounds on B 2 (E τ , ρ τ ), B 2 ( Ẽτ , ρτ ) and B 2 ( Êτ , ρτ ) pass to the limit by semicontinuity and allow to conclude that E, Ẽ and Ê are vector valued Radon measures absolutely continuous w.r.t. ρ. Hence there exist v t , ṽt , vt such that E = ρv, Ẽ = ρṽ and Ê = ρv.

Step 4. We now look at the equations satisfied by E, Ẽ and Ê. First we use ∂ t ρ τ + ∇ • E τ = 0, we pass to the limit as τ → 0, and we get

∂ t ρ + ∇ • E = 0.
Then, we use Ẽτ = -∇ρ τ + u t ρτ , we pass to the limit again as τ → 0, and we get Ẽ = -∇ρ + u t ρ.

To justify the above limit, the only delicate point is passing to the limit the term u t ρτ , since u is only L ∞ , and ρτ converges weakly as measures, and we are a priori only allowed to multiply it by continuous functions. Yet, we remark that by Corollary 3.3 we have that E(ρ τ t ) ≤ Cτ for all t ∈ [0, T ]. In particular, this provides, for each t, uniform integrability for ρτ t and turns the weak convergence as measures into weak convergence in L 1 . This allows to multiply by u t in the weak limit.

Finally, we look at Êτ . There exists a piecewise constant (in time) function p τ (defined as p τ k+1 on every interval ]kτ, (k + 1)τ ]) such that p τ ≥ 0, p τ (1 -ρτ ) = 0,

(3.6) T 0 Ω |∇p τ | 2 ( dx, dt) = T 0 Ω |∇p τ | 2 ρτ ( dx, dt) = T 0 Ω |v τ | 2 ρτ ( dx, dt) ≤ C
and Êτ = ∇p τ ρτ = ∇p τ . The bound (3.6) implies that p τ is uniformly bounded in L 2 (0, T ; H 1 (Ω)).

Since for every t we have |{p τ t = 0}| ≥ |{ρ τ t < 1}| ≥ |Ω| -1, we can use a suitable version of Poincaré's inequality, and get a uniform bound in L 2 ([0, T ]; L 2 (Ω)) = L 2 ([0, T ] × Ω). Hence there exists p ∈ L 2 ([0, T ] × Ω) such that p τ ⇀ p weakly in L 2 as τ → 0. In particular we have Ê = ∇p. Moreover it is clear that p ≥ 0 and by Lemma 3.5 we obtain p(1 -ρ) = 0 a.e. as well. Indeed, the assumptions of the Lemma are easily checked: we only need to estimate W 2 (ρ τ (a), ρτ (b)) for b > a, but we have

W 2 (ρ τ (a), ρτ (b)) = W 2 (ρ τ (k a τ ), ρ τ (k b τ )) ≤ C k b -k a , for k b τ ≤ b + τ and k a ≥ a.
Once we have Ê = ∇p with p(1 -ρ) = 0, p ∈ L 2 ([0, T ]; H 1 (Ω)) and ρ ∈ L ∞ , we can also write

Ê = ∇p = ρ∇p.
If we sum up our results, using E = Ẽ -Ê, we have

∂ t ρ -∆ρ + ∇ • (ρ(u -∇p)) = 0 together with p ≥ 0, ρ ≤ 1, p(1 -ρ) = 0 a.e. in [0, T ] × Ω.
As usual, this equation is satisfied in a weak sense, with no-flux boundary conditions.

Uniform Lip([0, T ]; W 1 ) and BV estimates

In this section we provide uniform estimates for the curves ρ τ , ρτ and ρτ of the following form: we prove uniform BV (in space) bounds on ρτ (which implies the same bound for ρτ ) and uniform Lipschitz bounds in time for the W 1 distance on ρ τ . This means a small improvement compared to the previous section in what concerns time regularity, as we have Lipschitz instead of AC 2 , even if we need to replace W 2 with W 1 . It is also important in what concerns space regularity. Indeed, from Lemma 3.2 one could deduce that the solution ρ of the FP equation (1.5) satifies √ ρ ∈ L 2 ([0, T ]; H 1 (Ω)) and, using ρ ≤ 1, also ρ ∈ L 2 ([0, T ]; H 1 (Ω)). Yet, this is just an integrable estimate in t, while the BV estimate of this section is uniform in the time variable.

Nevertheless there is a price to pay for this improvement: we have to assume higher regularity for the velocity field. These uniform-in-time W 1 -Lipschitz bounds are based both on BV estimates for the Fokker-Planck equation (see Lemma A.1 from Appendix A) and for the projection operator P K (see [START_REF] De Philippis | BV estimates in optimal transportation and applications[END_REF]). The assumption on u is essentially the following: we need to control the growth of the total variation of the solutions of the Fokker-Planck equation (3.1), and we need to iterate this bound along time steps.

We will discuss in the Appendix the different BV estimates on the Fokker-Planck equation that we were able to find. The desired estimate is true whenever u t C 1,1 (Ω) is uniformly bounded u t • n = 0 on ∂Ω. It seems to be an open problem to obtain similar estimate under the only assumption that u is Lipschitz continuous. Of course, we will also assume ρ 0 ∈ BV (Ω). Despite these extra regularity assumptions, we think these estimates have their own interest, exploiting some finer properties of the solutions of the Fokker-Planck equation and of the Wasserstein projection operator.

Before entering into the details of the estimates, we want to discuss why we concentrate on BV estimates (instead of Sobolev ones) and on W 1 (instead of W p , p > 1). The main reason is the role of the projection operator: indeed, even if ρ ∈ W 1,p (Ω), we do not have in general P K [ρ] ∈ W 1,p because the projection creates some jumps at the boundary of {P K [ρ] = 1}. This prevents from obtaining any W 1,p estimate for p > 1. On the other hand, [START_REF] De Philippis | BV estimates in optimal transportation and applications[END_REF] exactly proves a BV estimate on P K [ρ] and paves the way to BV bounds for our equation. Concerning the regularity in time, we observe that the velocity field in the Fokker-Planck equation contains a term in ∇ρ/ρ. Since the metric derivative in W p is given by the L p norm (w.r.t. ρ t ) of the velocity field, it is clear that estimates in W p for p > 1 would require spatial W 1,p estimates on the solution itself, which are impossible for p > 1 in this splitting scheme. We underline that this does not mean that uniform W 1,p are impossible for the solution of (1.5); it only means that they are not uniform along the approximation that we used in our Main Scheme to build such a solution.

The precise result that we prove is the following.

Theorem 4.1. Let us suppose that u t C 1,1 ≤ C and ρ 0 ∈ BV (Ω). Then using the notations from the Main scheme and Theorem 3.1 one has ρτ t BV ≤ C and W 1 (ρ τ k , ρ τ k+1 ) ≤ Cτ . As a consequence we also have

ρ ∈ Lip([0, T ]; W 1 ) ∩ L ∞ ([0, T ]; BV (Ω)).
To prove this theorem we need the following lemmas. Lemma 4.2. Suppose u t Lip ≤ C and u t • n = 0 on ∂Ω. Then for the solution ̺ of (A.1) with velocity field v = u we have the estimate

̺ t L ∞ ≤ ̺ 0 L ∞ e Ct , where C = ∇ • u t L ∞ .
Proof. Standard comparison theorems for parabolic equations allow to prove the results once we notice that f (t, x) := ̺ 0 L ∞ e Ct is a supersolution of the Fokker-Planck equation, i.e.

∂ t f t ≥ ∆f t -∇ • (f t u t ).

Indeed, in the above equation the Laplacian term vanishes as

f is constant in x, ∂ t f t = Cf t and ∇ • (f t u t ) = f t ∇ • u t + ∇f t • u t = f t ∇ • u t ≤ Cf t where C = ∇ • u t L ∞ .
From this inequality, and from ρ 0 ≤ f 0 , we deduce ρ t ≤ f t for all t.

We remark that the above lemma implies in particular that after every step in the Main scheme we have ρτ k+1 ≤ e τ c ≤ 1 + Cτ, where c := ∇ • u L ∞ . Let us now present the following lemma as well.

Corollary 4.3. Along the iterations of our Main scheme, for every k we have W 1 (ρ τ k+1 , ρ τ k+1 ) ≤ τ C for a constant C > 0 independent of τ .

Proof. With the saturation property of the projection (see Section 2.2 or [START_REF] De Philippis | BV estimates in optimal transportation and applications[END_REF]), we know that there exists a measurable set B ⊆ Ω such that ρ τ k+1 = ρτ k+1 ½ B + ½ Ω\B . On the other hand we know that

W 1 (ρ τ k+1 , ρ τ k+1 ) = sup f ∈Lip 1 (Ω), 0≤f ≤diam(Ω) Ω f (ρ τ k+1 -ρ τ k+1 ) dx = sup f ∈Lip 1 (Ω), 0≤f ≤diam(Ω) Ω\B f (ρ τ k+1 -1) dx ≤ τ C |Ω|diam(Ω).
We used the fact that the competitors f in the dual formula can be taken positive and bounded by the diameter of Ω, just by adding a suitable constant. This implies as well that C is depending on c, |Ω| and diam(Ω).

Proof of Theorem 4.1. First we take care of the BV estimate. Lemma A.1 in the Appendix guarantees, for t ∈]kτ, (k + 1)τ [, that we have T V (ρ τ t ) ≤ Cτ + e Cτ T V (ρ τ k ). Together with the BV bound on the projection that we presented in Section 2.2 (taken from [START_REF] De Philippis | BV estimates in optimal transportation and applications[END_REF]), this can be iterated, providing a uniform bound (depending on T V (ρ 0 ), T and sup t u t C 1,1 ) on ρτ t BV . Passing this estimate to the limit as τ → 0 we get ρ ∈ L ∞ ([0, T ]; BV (Ω)).

Then we estimate the behavior of the interpolation curve ρτ in terms of W 1 . We estimate

W 1 (ρ τ k , ρτ k+1 ) ≤ (k+1)τ kτ |(ρ τ t ) ′ | W 1 dt ≤ (k+1)τ kτ Ω |∇ρ τ t | ρτ t + |u t | ρτ t dx dt ≤ (k+1)τ kτ ρτ t BV dt + Cτ ≤ Cτ.
Hence, we obtain

W 1 (ρ τ k , ρ τ k+1 ) ≤ W 1 (ρ τ k , ρτ k+1 ) + W 1 (ρ τ k+1 , ρ τ k+1 ) ≤ τ C. This in particular means, for b > a, W 1 (ρ τ (a), ρτ (b)) ≤ C(b -a + τ ).
We can pass this relation to the limit, using that, for every t, we have ρτ t → ρ t in W 2 (Ω) (and hence also in

W 1 (Ω), since W 1 ≤ W 2 ), we get W 1 (ρ(a), ρ(b)) ≤ C(b -a),
which means that ρ is Lipschitz continuous in W 1 (Ω).

Variations on a theme: some reformulations of the Main scheme

In this section we propose some alternative approaches to study the problem (1.5). The general idea is to discretize in time, and give a way to produce a measure ρ τ k+1 starting from ρ τ k . Observe that the interpolations that we proposed in the previous sections ρ τ , ρτ and ρτ are only technical tools to state and prove a convergence result, and the most important point is exactly the definition of ρ τ k+1 .

The alternative approaches proposed here explore different ideas, more difficult to implement than the one that we presented in Section 3, and/or restricted to some particular cases (for instance when u is a gradient). They have their own modeling interest and this is the main reason justifying their sketchy presentation. 5.1. Variant 1: transport, diffusion then projection. We recall that the original splitting approach for the equation without diffusion ( [START_REF] Maury | Handling congestion in crowd motion modeling[END_REF][START_REF] Roudneff-Chupin | Modélisation macroscopique de mouvements de foule[END_REF]) exhibited an important difference compared to what we did in Section 3. Indeed, in the first phase of each time step (i.e. before the projection) the particles follow the vector field u and ρτ k+1 was not defined as the solution of a continuity equation with advection velocity given by u t , but as the image of ρ τ k via a straight-line transport: ρτ k+1 := (id + τ u kτ ) # ρ τ k . One can wonder whether it is possible to follow a similar approach here. A possible way to proceed is the following: take a random variable X distributed according to ρ τ k , and define ρτ k+1 as the law of X + τ u kτ (X) + B τ , where B is a Brownian motion, independent of X. This exactly means that every particle moves starting from its initial position X, following a displacement ruled by u, but adding a stochastic effect in the form of the value at time τ of a Brownian motion. We can check that this means

ρτ k+1 := η τ * ((id + τ u kτ ) # ρ τ k ) ,
where η τ is a Gaussian kernel with zero-mean and variance τ , i.e. η τ (x) := 1 (4τ π) d/2 e -|x| 2 4τ .

Then we define ρ τ k+1 := P K [ρ k+1 ] . Despite the fact that this scheme is very natural and essentially not that different from the Main scheme, we have to be careful with the analysis. First we have to quantify somehow the distance W p (ρ τ k , ρτ k+1 ) for some p ≥ 1 and show that this is of order τ in some sense. Second, we need to be careful when performing the convolution with the heat kernel (or adding the Brownian motion, which is the same): this requires either to work in the whole space (which was not our framework) or in a periodic setting (Ω = T d , the flat torus, which is qutie restrictive). Otherwise, the "explicit" convolution step should be replaced with some other construction, such as following the Heat equation (with Neumann boundary conditions) for a time τ . But this brings back to a situation very similar to the Main scheme, with the additional difficulty that we do not really have estimates on (id + τ u kτ ) # ρ τ k .

5.2. Variant 2: gradient flow techniques for gradient velocity fields. In this section we assume that the velocity field of the population is given by the opposite of the gradient of a function, u t = -∇V t a typical example is given when we take for V the distance function to the exit (see the discussions in [START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF] about this type of question). We start from the case where V does not depend on time, and we suppose V ∈ W 1,1 (Ω). In this particular case -beside the splitting approachthe problem has a variational structure, hence it is possible to show the existence by the means of gradient flows in Wasserstein spaces.

Since the celebrated paper of Jordan, Kinderlehrer and Otto ( [START_REF] Jordan | The variational formulation of the Fokker-Plack equation[END_REF]) we know that the solutions of the Fokker-Planck equation (with a gradient vector field) can be obtained with the help of the gradient flow of a perturbed entropy functional with respect to the Wasserstein distance W 2 . This formulation of the JKO scheme was also used in [START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF] for the first order model with density constraints. It is easy to combine the JKO scheme with density constraints to study the second order/diffusive model. As a slight modification of the model from [START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF], we can consider the following discrete implicit Euler (or JKO) scheme. As usual, we fix a time step τ > 0, ρ τ 0 = ρ 0 and for all k ∈ {1, 2, . . . , ⌊N/τ ⌋} we just need to define ρ τ k+1 . We take (5.1)

ρ τ k+1 = argmin ρ∈P(Ω) Ω V (x)ρ(x) dx + E(ρ) + I K (ρ) + 1 2τ W 2 2 (ρ, ρ τ k ) ,
where I K is the indicator function of K, which is

I K (x) := 0, if x ∈ K, +∞, otherwise.
The usual techniques from [START_REF] Jordan | The variational formulation of the Fokker-Plack equation[END_REF][START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF] can be used to identify that System (1.5) is the gradient flow of the functional ρ → J(ρ) := Ω V (x)ρ(x) dx + E(ρ) + I K (ρ) and that the above discrete scheme converges (up to a subsequence) to a solution of (1.5), thus proving existence. The key estimate for compactness is

1 2τ W 2 2 (ρ τ k+1 , ρ τ k ) ≤ J(ρ τ k ) -J(ρ τ k+1
), which can be summed up (as on the r.h.s. we have a telescopic series), thus obtaining the same bounds on B 2 that we used in Section 3.

Note that whenever D 2 V ≥ λI, the functional ρ → Ω V (x)ρ(x) dx+E(ρ)+I K (ρ) is λ-geodesically convex. This allows to use the theory in [START_REF] Ambrosio | Gradient flows in metric spaces and in the space of probability measure[END_REF] to prove not only existence, but also uniqueness for this equation, and even stability (contractivity or exponential growth on the distance between two solutions) in W 2 . Yet, we underline that the techniques of [START_REF] Marino | Uniqueness issues for evolution equations with density constraints[END_REF] also give the same result. Indeed, [START_REF] Marino | Uniqueness issues for evolution equations with density constraints[END_REF] contains two parts. In the first part, the equation with density constaints for a given velocity field u is studied, under the assumption that -u has some monotonicity properties: (-u t (x) + u t (y)) • (x -y) ≥ λ|x -y| 2 (which is the case for the gradients of λ-convex functions). In this case standard Grönwall estimates on the W 2 distance between two solutions are proved, and it is not difficult to add diffusion to that result (as the Heat kernel is already contractant in W 2 ). In the second part, via different techniques (mainly using the adjoint equation, and proving somehow L 1 contractivity), the uniqueness result is provided for arbitrary L ∞ vector fields u, but with the crucial help of the diffusion term in the equation.

It is also possible to study a variant where V depends on time. We assume for simplicity that V ∈ Lip([0, T ] × Ω) (this is a simplification; less regularity in space, such as W 1,1 , could be sufficient). In this case we define

J t (ρ) := Ω V t (x)ρ(x) dx + E(ρ) + I K (ρ) and (5.2) ρ τ k+1 = argmin ρ∈P(Ω) J kτ (ρ) + 1 2τ W 2 2 (ρ, ρ τ k ) ,
The analysis proceeds similarly, with the only exception that the we get

1 2τ W 2 2 (ρ τ k+1 , ρ τ k ) ≤ J kτ (ρ τ k ) -J kτ (ρ τ k+1 ),
which is no more a a telescopic series. Yet, we have J kτ (ρ τ k+1 ) ≥ J (k+1)τ (ρ τ k+1 ) + Lip(V )τ , and we can go on with a telescopic sum plus a remainder of the order of τ . In the case where u t is the opposite of the gradient of a λ-convex function V t , one could consider approximation by functions which are piecewise constant in time and use the standard theory of gradient flows.

Let us remark here that the recent paper [START_REF] Alexander | Quasi-static evolution and congested crowd transport[END_REF] gives another approach to deal with first order crowd motion models as limit of nonlinear-diffusion equations with gradient drift. This approach could be plausible also in the case when we add a simple diffusion term in the models studied in [START_REF] Alexander | Quasi-static evolution and congested crowd transport[END_REF]. 5.3. Variant 3: transport then gradient flow-like step with the penalized entropy functional. We present now a different scheme, which combines some of the previous approaches. It could formally provide a solution of the same equation, but presents some extra difficulties.

We define now ρτ k+1 := (id + τ u kτ ) # ρ τ k and with the help of this we define

ρ τ k+1 := argmin ρ∈K E(ρ) + 1 2τ W 2 2 (ρ, ρτ k+1 ).
In the last optimization problem we minimize a strictly convex and l.s.c. functionals, and hence we have existence and uniqueness of the solution. The formal reason for this scheme being adapted to the equation is that we perform a step of a JKO scheme in the spirit of [START_REF] Jordan | The variational formulation of the Fokker-Plack equation[END_REF] (without the density constraint) or of [START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF] (without the entropy term). This should let a term -∆ρ -∇ • (ρ∇p) appear in the evolution equation. The term ∇ • (ρu) is due to the first step (the definition of ρτ k+1 ). To explain a little bit more for the unexperienced reader, we consider the optimality conditions for the above minimization problem. Following [START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF], we can say that ρ ∈ K is optimal if and only if there exists a constant ℓ ∈ R and a Kantorovich potential ϕ for the transport from ρ to ρ τ k such that

ρ =      1 on ln ρ + ϕ τ < ℓ, 0 on ln ρ + ϕ τ > ℓ, ∈ [0, 1] on ln ρ + ϕ τ = ℓ.
We then define p = (ℓ -ln ρ -ϕ τ ) + and we get p ∈ press(ρ). Moreover, ρ-a.e. ∇p = -∇ρ ρ -∇ϕ τ . We then use the fact that the optimal transport is of the form T = id-∇ϕ and obtain a situation as is sketched in Figure 2.

• ρ τ k i d + τ u k τ • ρτ k+1 i d + τ ( ∇ p + ∇ ρ ρ ) • ρ τ k+1 id -τ (u (k+1)τ -∇p-∇ρ ρ )+o(τ )
Figure 2. One time step Notice that (id + τ u kτ ) -1 • (id + τ (∇p + ∇ρ/ρ)) = id -τ (u (k+1)τ -∇p -∇ρ/ρ) + o(τ ) provided u is regular enough. Formally we can pass to the limit τ → 0 and have

∂ t ρ -∆ρ + ∇ • (ρ(u -∇p)) = 0.
Yet, this turns out to be quite naïve, because we cannot get proper estimates on W 2 (ρ τ k , ρ τ k+1 ). Indeed, this is mainly due to the hybrid nature of the scheme, i.e. a gradient flow for the diffusion and the projection part on one hand and a free transport on the other hand. The typical estimate in the JKO scheme comes from the fact that one can bound W 2 (ρ τ k , ρ τ k+1 ) 2 /τ with the opposite of the increment of the energy, and that this gives rise to a telescopic sum. Yet, this is not the case whenever the base point for a new time step is not equal to the previous minimizer. Moreover, the main difficulty here is the fact that the energy we consider implicitly takes the value +∞, due to the constraint ρ ∈ K, and hence no estimate is possible whenever ρτ k+1 / ∈ K. As a possible way to overcome this difficulty, one could approximate the discontinuous functional I K with some finite energies of the same nature (for instance power-like entropies, even if the best choice would be an energy which is Lipschitz for the distance W 2 ). These kinds of difficulties are matter of current study, in particular for mixed systems and/or multiple populations.

Appendix A. BV -type estimates for the Fokker-Planck equation

Here we present some Total Variation (T V ) decay results (in time) for the solutions of the Fokker-Planck equation. Some are very easy, some trickier. The goal is to look at those estimates which can be easily iterated in time and combined with the decay of the T V via the projection operator, as we did in Section 4.

Let us take a vector field v : [0, +∞[×Ω → R d (we will choose later which regularity we need) and consider in Ω the problem

(A.1)    ∂ t ρ t -∆ρ t + ∇ • (ρ t v t ) = 0, in ]0, +∞[×Ω, ρ t (∇ρ t -v t ) • n = 0, on [0, +∞[×∂Ω, ρ(0, •) = ρ 0 ,
in Ω, for ρ 0 ∈ BV (Ω) ∩ P(Ω).

Lemma A.1. Suppose v t C 1,1 ≤ C for all t ∈ [0, +∞[. Suppose that either Ω = T d , or that Ω is convex and v • n = 0 on ∂Ω. Then, we have the following total variation decay estimate

(A.2) Ω |∇ρ t | dx ≤ C(t -s) + e C(t-s) Ω |∇ρ s | dx, ∀ 0 ≤ s ≤ t,
where C > 0 is a constant depending just on the C 1,1 norm of v.

Proof. First we remark that by the regularity of v the quantity v

L ∞ + Dv L ∞ + ∇(∇ • v) L ∞
is uniformly bounded. Let us drop now the dependence on t in our notation and calculate in coordinates

d dt Ω |∇ρ| dx = Ω ∇ρ |∇ρ| • ∇(∂ t ρ) dx = Ω ∇ρ |∇ρ| • ∇(∆ρ -∇ • (vρ)) dx = Ω j ρ j |∇ρ| i ρ iij -(∇ • (vρ)) j dx = - Ω i,j,k ρ 2 ij |∇ρ| - ρ j ρ k ρ ki ρ ij |∇ρ| 3 dx + B 1 - Ω j,i ρ j |∇ρ| v i ij ρ + v i i ρ j + v i j ρ i + v i ρ ij dx ≤ B 1 + C + C Ω |∇ρ| dx + Ω |∇ρ||∇ • v| dx + B 2 ≤ B 1 + B 2 + C + C Ω |∇ρ| dx.
Here the B i are the boundary terms, i.e.

B 1 := ∂Ω i,j ρ j n i ρ ij |∇ρ| dH d-1 and B 2 := - ∂Ω (v • n)|∇ρ| dH d-1 . The constant C > 0 only depends on v L ∞ + ∇ • v L ∞ + ∇(∇ • v) L ∞ .
We used as well the fact that -

Ω i,j,k ρ 2 ij |∇ρ| - ρ j ρ k ρ ki ρ ij |∇ρ| 3 dx ≤ 0.
Now, it is clear that in the case of the torus the boundary terms B 1 and B 2 do not exist, hence we conclude by Grönwall's lemma. In the case of the convex domain we have B 2 = 0 (because of the assumption v • n = 0) and B 1 ≤ 0 because of the next Lemma A.2.

Lemma A.2. Suppose that u : Ω → R d is a smooth vector field with u • n = 0 on ∂Ω, ρ is a smooth function with ∇ρ • n = 0 on ∂Ω, and that Ω ⊂ R d is a smooth convex set that we write as Ω = {h < 0} for a smooth convex function h with |∇h| = 1 on ∂Ω (so that n = ∇h on ∂Ω). Then we have, on the whole boundary ∂Ω, i,j

u i j ρ j n i = - i,j u i h ij ρ j .
In particular, we have i,j ρ ij ρ j n i ≤ 0.

Proof. The Neumann boundary assumption on u means u(γ(t)) • ∇h(γ(t)) = 0 for every curve γ valued in ∂Ω and for all t. Differentiating in t, we get i,j

u i j (γ(t))(γ ′ (t)) j h i (γ(t)) + i,j u i (γ(t))h ij (γ(t))(γ ′ (t)) j = 0.
Take a point x 0 ∈ ∂Ω and choose a curve γ with γ(t 0 ) = x 0 and γ ′ (t 0 ) = ∇ρ(x 0 ) (which is possible, since this vector is tangent to ∂Ω by assumption). This gives the first part of the statement. The second part, i.e.

i,j ρ ij ρ j n i ≤ 0, is obtained by taking u = ∇ρ and using that D 2 h(x 0 ) is a positive definite matrix.

Remark A.3. If we look attentively at the proof of Lemma A.1, we can see that we did not really exploit the regularizing effects of the diffusion term in the equation. This means that the regularity estimate that we provide are the same that we would have without diffusion: in this case, the density ρ t is obtained from the initial density as the image through the flow of v. Thus, the density depends on the determinant of the Jacobian of the flow, hence on the derivatives of v. It is normal that, if we want BV bounds on ρ t , we need assumptions on two derivatives of v.

We would like to prove some form of BV estimates under weaker regularity assumptions on v, trying to exploit the diffusion effects. In particular, we would like to treat the case where v is only C 0,1 . As we will see in the following lemma, this degenerates in some sense.

Lemma A.4. Suppose that Ω is either the torus or a smooth convex set Ω = {h < 0} parameterized as a level set of a smooth convex function h. Let v t : Ω → R d be a vector field for t ∈ [0, T ], Lipschitz and bounded in space, uniformly in time. In the case of a convex domain, suppose v • n = 0 on ∂Ω. Let H : R d → R be given by H(z) := ε 2 + |z| 2 . Now let ρ t (sufficiently smooth) be the solution of the Fokker-Planck equation with homogeneous Neumann boundary condition.

Then there exists a constant C > 0 (depending on v and Ω) such that

(A.3) Ω H(∇ρ t ) dx ≤ Ω H(∇ρ 0 ) dx + Cεt + C ε t 0 ρ s 2 L ∞ ds.
Proof. First let us discuss about some properties of H. It is smooth, its gradient is

∇H(z) = z H(z) and it satisfies ∇H(z) • z ≤ H(z), ∀z ∈ R d . Moreover its Hessian matrix is given by [H ij (z)] i,j∈{1,...,d} = δ ij H 2 (z) -z i z j H 3 (z) i,j∈{1,...,d} = 1 H(z) I d - 1 H 3 (z) z ⊗ z, ∀ z ∈ R d , where δ ij = 1, if i = j, 0, if j = j,
is the Kronecker symbol. Note that, from this computation, the matrix D 2 H ≥ 0 is bounded from above by 1 H , and hence by ε -1 . Moreover we introduce a Now we study each term separately and for the simplicity we drop the t subscripts in the followings. We start from the case of the torus, where there is no boundary term in the integration by parts. First look at the term (II a ). Since the matrix H jk is positive definite, we can apply a Young inequality for each index i and obtain

uniform constant C > 0 such that v 2 L ∞ |Ω| + ∇ • v L ∞ + Dv L ∞ ≤ C. Now to
(II a ) = Ω i,j,k H jk (∇ρ)ρ ki v i j ρ dx ≤ 1 2 Ω i,j,k H jk (∇ρ)ρ ki ρ ij dx + 1 2 Ω i,j,k H jk (∇ρ)v i j v i k ρ 2 dx ≤ 1 2 |(I)| + C ρ 2 L 2 D 2 H L ∞ .
The L 2 norm in the second term will be estimated by the L ∞ norm for the sake of simplicity (see Remark A.5 below).

For the term (II b ) we first make a point-wise computation i,j,k

H jk (∇ρ)ρ ki v i ρ j = 1 H 3 (∇ρ) i [D 2 i ρ • ε 2 I d + |∇ρ| 2 I d -∇ρ ⊗ ∇ρ • ∇ρ]v i = ε 2 H 3 (∇ρ) i v i D 2 i ρ • ∇ρ = -ε 2 i v i ∂ i 1 H(∇ρ) .
where D 2 i ρ denotes the i th row in the Hessian matrix of ρ and we used |∇ρ| 2 I d -∇ρ ⊗ ∇ρ •∇ρ = 0. Integrating by parts we obtain

(II b ) = ε 2 Ω (∇ • v) 1 H(∇ρ) dx ≤ Cε 2 1/H L ∞ ≤ Cε,
where we used H(z) ≥ ε.

Summing up all the terms we get and using D 2 H ≤ ε -1 we get

d dt Ω H(∇ρ t ) dx ≤ - 1 2 |(I)| + C ρ t 2 L ∞ D 2 H L ∞ + Cε ≤ Cε + C ρ t 2 L ∞ ε -1 ,
which proves the claim.

If we switch to the case of a smooth bounded convex domain Ω, we have to handle boundary terms. These terms are ∂Ω i,j H j (∇ρ)ρ ij n i -∂Ω i,j H j (∇ρ)ρv i j n i , where we ignored those terms involving n i v i (i.e., the integration by parts in (II b ), and the term H j (∇ρ)ρ j n i v i in the integration by parts of (II a )), since we already supposed v • n = 0. We use here Lemma A.2, which provides i,j

H j (∇ρ)ρ ij n i -ρH j (∇ρ)v i j n i = 1 H(∇ρ) i,j ρ j ρ ij n i -ρρ j v i j n i = -1 H(∇ρ) i,j ρ j h ij ρ i -ρρ j h ij v i .

If we use the fact that the matrix D 2 h is positive definite and a Young inequality, we get i,j ρ j h ij ρ i ≥ 0 and

ρ i,j |ρ j h ij v i | ≤ 1 2 i,j ρ j h ij ρ i + 1 2 i,j ρ 2 v j h ij v i , which implies 1 H(∇ρ) i,j ρ j ρ ij n i -ρρ j v i j n i ≤ ρ 2 H(∇ρ) D 2 h L ∞ |v| 2 ≤ C ρ 2 L ∞ ε .
This provides the desired estimate on the boundary term.

Remark A.5. In the above proof, we needed to use the L ∞ norm of ρ only in the boundary term. When there is no boundary term, the L 2 norm is enough, in order to handle the term (II a ). In both cases, the norm of ρ can be bounded in terms of the initial norm multiplied by e Ct , where C bounds the divergence of v. On the other hand, in the torus case, one only needs to suppose ρ 0 ∈ L 2 and in the convex case ρ 0 ∈ L ∞ . Both assumptions are satisfied in the applications to crowd motion with density constraints.

We have seen that the constants in the above inequality depend on ε and explode as ε → 0. This prevents us to obtain a clean estimate on the BV norm in this context, but at least proves that ρ 0 ∈ BV ⇒ ρ t ∈ BV for all t > 0 (to achieve this result, we just need to take ε = 1). Unfortunately, the quantity which is estimated is not the BV norm, but the integral Ω H(∇ρ). This is not enough for the purpose of the applications to Section 4, as it is unfortunately not true that the projection operator decreases the value of this other functional 2 .

If we stay interested to the value of the BV norm, we can provide the following estimate.

Lemma A.6. Under the assumptions of Lemma A.4, if we suppose ρ 0 ∈ BV (Ω) ∩ L ∞ (Ω), then, for t ≤ T , we have

(A.4) Ω |∇ρ t | dx ≤ Ω |∇ρ 0 | dx + C √ t,
where the constant C depends on v, on T and on ρ 0 L ∞ .

2 Here is a simple counter-example: consider µ = g(x) dx a BV density on [0, 2] ⊂ R, with g defined as follows.

Divide the interval [0, 2] into 2K intervals Ji of length 2r (with 2rK = 1); call ti the center of each interval Ji (i.e. ti = i2r + r, for i = 0, . . . , 2K -1) and set g(x) = L + r 2 -(x -ti) 2 on each Ji with i odd, and g(x) = 0 on Ji for i even, taking L = 1 -πr/4. It is not difficult to check that the projection of µ is equal to the indicator function of the union of all the intervals Ji with i odd, and that the value of H(∇ρ) has increased by K(2 -π/2)r = 1 -π/4, i.e. by a positive constant (see Figure 3). It is sufficient to choose, for fixed t, ε = √ t, in order to prove the claim.

Unfortunately, this √ t behavior is not suitable to be iterated, and the above estimate is useless for the sake of Section 4. The existence of an estimate (for v Lipschitz) of the form T V (ρ t ) ≤ T V (ρ 0 ) + Ct, or T V (ρ t ) ≤ T V (ρ 0 )e Ct , or even f (T V (ρ t )) ≤ f (T V (ρ 0 ))e Ct , for any increasing function f : R + → R + , seems to be an open question.

Lemma 3 . 2 .

 32 Let us consider a solution ̺ t of the Fokker-Planck equation on [0, T ] × Ω with the velocity field u satisfying (U) and with no-flux boundary conditions on [0, T ] × ∂Ω. Then for any time interval ]a, b[ we have the following estimate

Corollary 3 . 3 .

 33 4 ([a, b] × Ω), keeping the L ∞ bound). If we call ̺ k the corresponding solution of the Fokker Planck equation, it satisfies (3.3). This implies a uniform bound (w.r.t. k) for ̺ k in L 2 ([a, b]; H 1 (Ω)), and hence a uniform bound on ̺ k in L 2 ([a, b] × Ω). From these bounds and the uniqueness of the solution of the Fokker-Planck equation with L ∞ drift we deduce ̺ k → ̺. The semicontinuity of the left-hand side in (3.3) and of the entropy term at t = b, together with the convergence of the entropy at t = a and the convergence b a Ω |u k | 2 ̺ k dx dt → b a Ω |u| 2 ̺ dx dt (because we have a product of weak and strong convergence in L 2 ) allow to pass (3.3) to the limit. From the inequality (3.4) we deduce that

for a C > 0

 0 independent of τ . Here we denoted by |(ρ τ t ) ′ | W 2 the metric derivative of the curve ρ τ at t in W 2 . In particular, we have a uniform Hölder bound on ρ τ : W 2 (ρ τ (a), ρ τ (b)) ≤ C √ b -a for every b > a. (iv) E τ , Ẽτ , Êτ are uniformly bounded sequences in M([0, T ] × Ω) d . Proof. (i) First by the triangle inequality and by the fact that ρ τ k+1

(

  ) • ∇∆ρ dx = Ω j,i H j (∇ρ)ρ jii dx = -Ω j,i,k H kj (∇ρ)ρ ik ρ ji dx (II) = -Ω ∇H(∇ρ) • ∇(∇ • (vρ)) dx = -Ω i,j H j (∇ρ)(v i ρ) ij dx = Ω i,j,k H jk (∇ρ)ρ ki v i j ρ dx + Ω i,j,kH jk (∇ρ)ρ ki v i ρ j dx =: (II a ) + (II b ).

Figure 3 .

 3 Figure 3. The counter-example to the decay of

  is (again) the solution of the Fokker-Planck equation (3.1) on the time interval [0, τ ] with initial datum ρ τ k . Here we do not double its speed. We define the curve of velocities ): 1) the first one is continuous in time for the weak-* convergence, while the second and third ones are not;

		ṽτ t := -	∇̺ t-kτ ̺ t-kτ	+ u t , if t ∈ [kτ, (k + 1)τ [ ,
	and we build the curve of momenta by Ẽτ t := ρτ t	ṽτ t .
	Third interpolation. For each τ , we also define piecewise constant curves,
		ρτ t := ρ τ k+1 , t := ∇p τ vτ k+1 ,	if t ∈ [kτ, (k + 1)τ [, if t ∈ [kτ, (k + 1)τ [,
	and Êτ t := ρτ t	vτ t . We remark that p τ k+1 (1 -ρ τ k+1 ) = 0, hence the curve of momenta is just
		Êτ
				τ [ ,
	and, finally, we simply define the curve of momenta as E τ t := ρ τ t v τ t .
	Second interpolation. We define another interpolation as follows. Set
		ρτ

t := ̺ t-kτ , if t ∈ [kτ, (k + 1)τ [, where ̺ t t := ∇p τ k+1 , if t ∈ [kτ, (k + 1)τ [. Mind the differences in the construction of ρ τ t , ρτ t and ρτ t (hence in the construction of v τ t , ṽτ t and vτ t and E τ t , Ẽτ t and Êτ t

  show the estimate of this lemma we calculate the quantity d dt H(∇ρ t ) dx.

	d dt Ω	H(∇ρ t ) dx =

Ω ∇H(∇ρ t ) • ∂ t ∇ρ t dx = Ω ∇H(∇ρ t ) • ∇(∆ρ t -∇ • (v t ρ t )) dx = Ω ∇H(∇ρ t ) • ∇∆ρ t dx -Ω ∇H(∇ρ t ) • ∇(∇ • (v t ρ t )) dx =: (I) + (II)

Acknowledgements. The authors warmly acknowledge the support of the ANR project ISO-TACE (ANR-12-MONU-0013) and of the iCODE project "strategic crowds", funded by the IDEX Université Paris-Saclay. They also acknowledge the warm hospitality of the Fields Institute of Toronto, where a large part of the work was accomplished during the Thematic Program on Variational Problems in Physics, Economics and Geometry in Fall 2014. Last but not least, the two anonymous referees who read very carefully the paper and suggested many useful improvements are warmly thanked as well.