Christophe Tournassat
email: c.tournassat@brgm.

Stéphane Gaboreau

Jean-Charles Robinet

Ian C Bourg

Carl I Steefel

Impact of microstructure on anion exclusion in compacted clay media

Keywords: clay, illite, montmorillonite, anion exclusion, Poisson-Boltzmann, microstructure, molecular dynamics


Pierre Buyssens

Isabelle Gardin

Su Ruan

Eikonal based region growing for superpixels generation : Application to semi-supervised real time organ segmentation in CT images

Keywords: Superpixels, Eikonal equation, semi-supervised segmentation

We describe a semi-supervised organ segmentation method for Computed Tomography images. In a first step, a dense oversegmentation of the image is created with an Eikonal-based algorithm. The proposed superpixel algorithm ourperforms state-of-the-art algorithms on classical metrics. In a second step, the semi-supervised segmentation is performed on the underlying Region Adjacency Graph created from the oversegmentation. As the complexity is greatly reduced, the organ segmentation can be performed in real-time.

Introduction

Fast and accurate computer aided organ segmentation is more and more required in many medical applications such automatic detection of liver cancer [START_REF] Shimizu | Preliminary report of cad system competition for liver cancer extraction from 3d ct images and fusion of the cads[END_REF], measurements of liver volume [START_REF] Nakayama | Automated hepatic volumetry for living related liver transplantation at multisection ct1[END_REF], helpful both for surgical planning prior to hepatic resection [START_REF] Fan | Safety of donors in live donor liver transplantation using right lobe grafts[END_REF], and for developing effective radiation treatments programs [START_REF] Kobashi | Knowledge-based organ identification from ct images[END_REF]. More precisely our work is part of a Selective Internal Radiation Therapy (SIRT) project [START_REF] Dieudonné | Clinical feasibility of fast 3-dimensional dosimetry of the liver for treatment planning of hepatocellular carcinoma with 90y-microspheres[END_REF] whose goal is to compute a 3dimensional dosimetry of the liver for treatment planning.

To this end the volume of the liver has to be known to deliver the right dose to the patient. Lung volume has also to be known to compute the Lung Shunt Fraction (LSF) in order to avoid irradiation of the lungs during the treatment [START_REF] Ahmadzadehfar | Radioembolization of liver tumors with yttrium-90 microspheres[END_REF]. In this work we propose a semi-supervised organ segmentation scheme that gives the control of the segmentation to the radiologist. A major challenge of such a scheme in medical applications is the volume of data that has to be processed, which is often incompatible with real-time processings.

To tackle this problem, we propose a twostep scheme that, first simplifies Computed Tomography images by computing an oversegmentation of the volumes, and second performs the segmentation on the graph created from this oversegmentation. The oversegmentation proposed in this paper belongs to the superpixel algorithms that aim at grouping perceptually and adjacent pixels into meaningful regions, the so-called superpixels.

Numerous superpixel algorithms have been proposed in the literature, they can roughly be divided into Graph-based methods [START_REF] Shi | Normalized cuts and image segmentation. Pattern Analysis and Machine Intelligence[END_REF][START_REF] Veksler | Superpixels and supervoxels in an energy optimization framework[END_REF] and gradient (or color) based ones [START_REF] Levinshtein | Turbopixels: Fast superpixels using geometric flows[END_REF][START_REF] Vedaldi | Quick shift and kernel methods for mode seeking[END_REF][START_REF] Achanta | Slic superpixels compared to state-of-the-art superpixel methods[END_REF]. More details can be found on these methods in Section 3 and a review of these in [START_REF] Achanta | Slic superpixels compared to state-of-the-art superpixel methods[END_REF]. Other works make use of similar oversegmentations such [START_REF] Andrés | Energy partitions and image segmentation[END_REF][START_REF] Arbeláez | A metric approach to vector-valued image segmentation[END_REF] based on an energy that favors partitions of similar pixels, or [START_REF] Ng | Medical image segmentation using k-means clustering and improved watershed algorithm[END_REF] based on the watershed transformation. These last methods however do not belong to the superpixel algorithm class since they do not control the amount of regions nor their compactness.

Eikonal-based Region Growing Clustering algorithm 1 (ERGC) is a new superpixel algorithm that formulates the clustering of pixels as a solution of an Eikonal equation :

∇U(x) = F(x) ∀x ∈ I U(x) = 0 ∀x ∈ Γ (1)
where I is the image domain, F a positive function, Γ the set of initial seeds, and U(x) the traveling time or geodesic distance of x from source Γ.

The main difference from other eikonalbased algorithms comes from the computation of F(x) which is based on a distance between pixel and region features.

Efficiently numericaly solved with the Fast Marching method [START_REF] Sethian | Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision[END_REF], ERGC leads to the creation of homogeneous superpixels. An additional spatial constraint can also be added to F to enforce superpixel compactness.

The rest of the paper is organized as follows : Section 2 details ERGC algorithm. Its performance is evaluated and compared to state-ofthe-art superpixel algorithms in Section 3. The application of ERGC as a preprocessing step of the graph-based semi-supervised segmentation scheme is detailed in Section 4. Section 5 concludes the paper. 1 Source code and executable of ERGC can be found at https://sites.google.com/site/ pierrebuyssens/ergc 2. Superpixels method

Notations

In the following several notations are adopted to simplify the reading of the paper. A pixel/voxel is noted p and consists of a coordinate triplet (x p , y p , z p ) and a color vector C p . For 2D images, the z coordinate is equal to zero. A superpixel/supervoxel R i consists of a seed pixel/voxel s i , a mean color vector C i and the number N i of pixels/voxels that belong to R i .

Note that these formulations are independant of the number of channels of the images. For gray level images, C p reduces to the gray level of the pixel p.

Algorithm

The initialization of the algorithm consists of sampling the number of desired seeds pixels/voxels K on a regular grid. Without any prior, the sampling step for 2D images is set to S = √ N/K with N the number of pixels/voxels of the image (similar to [START_REF] Achanta | Slic superpixels compared to state-of-the-art superpixel methods[END_REF] and [START_REF] Levinshtein | Turbopixels: Fast superpixels using geometric flows[END_REF]). These initial seeds are then perturbed and moved to the minima of the gradient map in a V 8 neighborhood (V 26 for 3D images). The mean color vector of each superpixel/supervoxel C i is then initialized with C S i and its neighbors in V 4 (V 6 for 3D images).

To compute the numerical solution of equation 1, we adopt the Fast Marching method introduced by Sethian in [START_REF] Sethian | Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision[END_REF]. It consists of computing locally the geodesic distances for an evolving front with a priority queue (see [START_REF] Peyré | Geodesic methods in computer vision and graphics[END_REF] for a detailed description of the fast marching method). As a front can only evolve with a positive speed function, it passes through a pixel/voxel only once, and the algorithm terminates when all pixels/voxels have been visited.

Contrary to many front-based methods that rely on a fixed potential map (often the magnitude of the gradient of the image : F = ∇I ), the right term of equation 1 evolves during the diffusion process and is not fixed. Given a pixel p, the proposed function F (right term of equation 1) is computed according to an adjoining region R i as :

F c ( p, R i ) = C p -C i 2 2
(2)

This formalization favors the creation of superpixels/supervoxels that are more homogeneous than those created by classical gradientbased methods. Figure 1 compares the behaviors of both approaches on a toy problem and on a natural scene.

Within the fast marching algorithm, a superpixel/supervoxel R i has to be updated every time a pixel p is integrated to R i . This update is easily computed with :

C i ← C i × N i + C p N i + 1 N i ← N i + 1
A similar potential function has been proposed in [START_REF] Li | Fast 3d brain segmentation using dual-front active contours with optional user-interaction[END_REF]. Dedicated to the 3D brain segmentation, the potential is computed according to a histogram analysis prior to the evolution of active contours. Contrary to our approach, the features (i.e.: mean and standard deviation of the color) of the regions that initially define the active contours do not evolve. Moreover the analysis of the histogram of the brain allows the detection of 3 well defined peaks (corresponding to the gray matter, the white matter and the cerebrospinal fluid), which is not the case for many kind of images, especially with CT where different organs may be composed of pixels with similar or equal gray values. Our framework is far more general since it can deal with any kind of multispectral images.

In this primal form, ERGC superpixels may be very irregular in shape and size. However an additional spatial constraint can be added to the right member of equation 1 to enforce superpixel compactness.

Superpixels compactness

According to the application, it may be desirable to have compact superpixels. Following the idea introduced in [1], a spatial constraint term can be added to the computation of F. This spatial constraint penalizes the geodesic distances of a pixel p that is spatially far from a seed s i via the L 2 norm ps i 2 .

Since the spatial and color proximities are not defined on the same range (especially with medical images), a compactness factor m is used to weigh the relative importance between the two metrics.

With this spatial constraint, the second term of equation 1 becomes :

F( p, R i ) = F c ( p, R i ) 2 + F s (p, s i ) 2 (3)
with

F s (p, s i ) = p -s i 2 2 S × m (4)
where S in the initial sampling step. Setting m = 0 leads to 2 which means that no spatial constraint is considered.

F( p, R i ) = F c ( p, R i )
Figure 2 shows the resulting superpixels for a CT image with varying the number of superpixels and the compactness factor.

3D extension

ERGC naturally extends to 3D images with minor modifications of the algorithm. Figure 3 shows a supervoxel oversegmentation of a 3D volume from the MICCAI-2007 Grand Challenge dataset. For display purposes only the interior part of the body is shown.

Evaluation of ERGC

In this section, we evaluate ERGC algorithm on a natural image dataset and on a medical image dataset, and compare it to the superpixel methods published in the litterature.

The classical metrics to compare oversegmentation results are used in this evaluation : • The boundary recall (BR) indicates how well superpixels adhere to true boundaries of the objects in the image. It measures the fraction of ground truth edges that is also present in superpixel segmentation within a distance t. In our experiments, t is fixed to 2 (as in [START_REF] Achanta | Slic superpixels compared to state-of-the-art superpixel methods[END_REF][START_REF] Levinshtein | Turbopixels: Fast superpixels using geometric flows[END_REF] to cite few).

• The under-segmentation (UE) error measures the fraction of bleeding caused by superpixels that overlap a given ground truth segment. For a given ground truth segmentation into segments g 1 , . . . , g k and a superpixel segmentation R 1 , . . . , R N , the undersegmentation error for a segment g i is quantified with

{R j |R j ∩g i >B} Area(R j ) -Area(g i ) Area(g i )
where B is set to 5 percent of Area(R j ) (as in [START_REF] Achanta | Slic superpixels compared to state-of-the-art superpixel methods[END_REF]) to account for ambiguities of the ground truth segments. Finally this fraction is simply averaged across all ground truth segments and all images.

• The Achievable Segmentation Accuracy (ASA) is a performance upperbound measure. It gives the highest performance when taking superpixels as units for object segmentation. Each superpixel is labeled with the label of the ground truth segment with the largest overlap. The fraction of correctly labeled pixels is the achievable accuracy :

AS A(s) = 100 × k max i |s k ∩ g i | i g i •
As superpixels are often used as a preprocessing step, the processing time is an important factor such that it does not slow down the image processing pipeline too much.

Performances are evaluated and compared to five other algorithms :

• The Normalized Cuts approach (NC05)

proposed in [START_REF] Cour | Spectral segmentation with multiscale graph decomposition[END_REF] partition a graph defined on pixels by minimizing a global cost function with an explicit compactness constraint.

• Superpixels method proposed in [START_REF] Veksler | Superpixels and supervoxels in an energy optimization framework[END_REF] formulates the superpixel partitioning problem in an energy minimization framework, and optimize it with graph cuts.

Two versions of the algorithm are proposed differing in compact (GCa10) or constant intensity (GCb10) resulting superpixels.

• The turbopixels method (TP09) proposed in [START_REF] Levinshtein | Turbopixels: Fast superpixels using geometric flows[END_REF] dilates a set of initial seeds using a level-set geometric flow to adapt to local image structure.

• Quick shift method (QS08) proposed in [START_REF] Vedaldi | Quick shift and kernel methods for mode seeking[END_REF] exploits kernels to extend mean shift [START_REF] Comaniciu | Mean shift: A robust approach toward feature space analysis[END_REF] using a medoid shift procedure.

• SLIC method proposed in [START_REF] Achanta | Slic superpixels compared to state-of-the-art superpixel methods[END_REF] adapts kmeans in the image plane to exhibit superpixels. The addition of a spatial constraint term produces regular regions that adhere well to image boundaries.

Berkeley dataset

The Berkeley image dataset contains 300 2D images (in the RGB colorspace) of size 481 × 312 (or 321 × 481) and approximately 10 ground truth manual segmentations for each image.

Table 1 exhibits results on standard metrics for the considered algorithms. Results have been obtained with approximately 500 superpixels per image and have been computed for the 3000 ground truth images. ERGC outperforms all other superpixels algorithms on the standard metrics, while being faster. Moreover, Student's t-test on these results exhibits high values (t > 35), which validates the statistical significance of the results.

Complexity of ERGC relies essentialy on the complexity of the fast marching algorithm, which is roughly in O(n log n) with an appropriate heap for sorting the pixels according to their geodesic distance. Despite this theorical complexity, ERGC is very fast in practice, and the algorithm is nearly linear in time. Note that these results have been obtained without spatial constraints m = 0. Adding this constraint degrades slightly ERGC performances, but produces more visually pleasing superpixels.

MICCAI-2007 Grand Challenge dataset

The MICCAI-2007 Grand Challenge dataset [START_REF] Heimann | Comparison and evaluation of methods for liver segmentation from ct datasets[END_REF] consists of 20 CT 3D images (volumes) used for training and 10 for tests. In this evaluation we focus on the training part (and its ground truth) and especially the slices which contain some liver parts. The sub-dataset consists then of 2750 2D gray images (of size 512 × 512), for a total amount of more than 1.5 million liver boundary pixels. Table 2 summarizes the results of ERGC and the current state-of-the-art SLIC algorithm. Since CT images are grayscale and for fair comparisons purposes, results of SLIC and ERGC of table 2 have been obtained by simulating RGB images (i.e. 3 channels) for a fair comparison. ERGC g is more realistic and works directly on monochannel images. ERGC g then gives the same results as ERGC while being faster.

As shown in table 1 and 2, ERGC and SLIC give the lowest Undersegmentation Er- The Achievable Segmentation Accuracy result is the more important in our application since the semi-supervised segmentation proposed below considers superpixels/supervoxels as units.

Semi-supervised organ segmentation

In order to simplify the images and to reduce the complexity of the semi-supervised segmentation algorithm, ERGC is performed on CT images. The segmentation is then processed on the underlying Region Adjacency Graph (RAG) created from this oversegmentation.

Feature of each vertex is simply the mean color of the underlying superpixel/supervoxel. An edge connecting two vertices u and v is weighted with the function

w(u, v) = C u -C v 2
Figure 4 shows a 512 × 512 slice extracted from a CT volume of the MICCAI-2007 Grand Challenge dataset, the oversegmentation performed by ERGC containing approximately 1000 superpixels, and the RAG created from this oversegmentation. In this example, no spatial constraint has been added to F, that is m = 0.

Once the graph is created, a user can manually label a superpixel and then its corresponding vertex in the graph. Each time a superpixel/supervoxel is labeled, a label diffusion on graph is performed using the Dijkstra algorithm [START_REF] Wybe | A short introduction to the art of programming[END_REF]. Since the graph contains very few vertices (compared to the number of pixels of the initial image), the diffusion is very fast. The labeled graph then gives a label to each superpixel/supervoxel, which leads to the segmentation.

In our tests, the Dijkstra diffusion on graph that contains 1000 vertices processes in approximately 0.01 seconds on a standard laptop. This processing time allows a real-time segmentation by the radiologist, and a bad labeling can be corrected on the fly. Figure 5 shows the semi-supervised segmentation of the liver on a 2D slice with few manual labels.

Processing 3D CT images operates in the same manner. Supervoxels are computed from a CT volume and the underlying RAG is created. However, instead of fixing the number of desired K, the initial seeds of ERGC are sampled according to a fix sampling step S . In our experiments on the MICCAI-2007 Grand Challenge dataset, S is fixed to 10 pixels along x and y axis, and 5 pixels along z axis, to take into account the inter-slices distance which is higher than inter-pixels distance in the transverse plane. The second advantage of fixing S along each axis instead of K is to deal with images of various depth.

Figure 6 shows transverse, sagittal and coronal views of a CT volume containing 183 slices of size 512 × 512. About 96000 ERGC supervoxels have been computed with the sampling steps detailed above. Fixing a high value for the spatial constraint m gives more compact supervoxels, hence simplifying the semi-supervised segmentation. In our experiments, m is fixed to the product of sampling steps which is 500.

Since the graph contains more vertices than in the 2D example above, the Dijkstra diffusion processing time is quite longer especially after the firsts supervoxels seeds have been placed. For the underlying RAG of figure 6, the first Dijkstra iteration processes in about 0.3 seconds but the following quickly decreases to 0.01 seconds or less, since distances for only a few number of vertices have to be updated.

Placing seeds into transverse, sagittal and coronal planes is natural for a radiologist. The whole process acts in the same manner as in 2D except that seeds can be placed in any of the three planes. Result segmentation is then updated in the three planes and the radiologist can add/remove seeds to refine/correct the segmentation. Figure 6 shows a segmented volume with liver, lungs and background belonging to different classes (i.e. having different labels). Only a few supervoxel seeds are needed to obtain such a result. Note that edges of lungs and liver are smoothed with gaussian kernels of different strengths.

Conclusion

In this paper, we proposed a method based on the Eikonal equation that quickly creates accurate superpixels/supervoxels. Outperforming other oversegmentation algorithms on classical metrics, ERGC is used as a first step in a graph-based semi-supervised image segmentation. Applied to Computed Tomography images, it allows a real-time segmentation by processing the underlying graph, and gives the control of the segmentation to the radiologist.

Further works will take into account color and shape prior to ease the organs segmentation.

Figure 1 :

 1 Figure 1: First row : initial image. The back and white circles depict the initial seeds. Second row : the result segmentation with the proposed F (left) and the gradient-based potential function (right). Third row : geodesic distances map U in fake color. Some isocontours are shown in white. Last row : segmentation of a natural scene. Seeds are depicted in black.

Figure 3 :

 3 Figure 3: Approximately 1000 ERGC supervoxels for a Computed Tomography volume of a human body. Only the interior part of the body is shown for visualization convenience.

(a) m = 0 (b) m = 10 (c) m = 100 Figure 2 :

 0101002 Figure 2: CT image with 100 superpixels in the upper left of each image and 900 in the lower right for different values of the compactness factor.

Figure 4 :

 4 Figure 4: (a) Initial slice, (b) approximately 1000 ERGC superpixels without spatial constraint (m = 0), and (c) the underlying RAG.

Figure 5 :

 5 Figure 5: Example of a semi-supervised segmentation of the liver. Two superpixels have been labeled as liver in dark green (a), (b,c,d,e,f) annotation of superpixels as background in dark red. The liver contour is computed on the fly and appears in green. For clarity purposes, only the liver is segmented.

Figure 6 :

 6 Figure 6: Segmentation of lungs (blue) and liver (green) in transverse (a), sagittal (b), and coronal (c) planes. Supervoxels seeds for the lungs, the liver and the background are colored dark blue, dark green and dark red respectively.

Table 2 :

 2 Boundary Recall, Undersegmentation Error, Achievable Segmentation Accuracy and processing time (in seconds) comparisons between SLIC, ERGC and ERGC g on the training part of the MICCAI-2007

	Grand Challenge dataset. Processing times are given
	for one slice.
		SLIC	ERGC ERGC g
	BR
	mean	0.8736 0.9464 0.9464
	std	0.0317 0.0209 0.0209
	UE	0.15	0.15	0.15
	ASA	98.03	99.12	99.12
	Processing time	1.56	0.57	0.37
	ror. However ERGC gives the better Bound-
	ary Recall and Achievable Segmentation Ac-
	curacy results upon the considered algorithms,
	and especially SLIC.

Table 1 :

 1 Comparison of ERGC to several superpixels algorithms for approximately 500 superpixels on different metrics. Values of algorithms NC05, GCa10, TP09, and QS08 have been extracted from[START_REF] Achanta | Slic superpixels compared to state-of-the-art superpixel methods[END_REF] courtesy of the authors.

			Graph-based			Gradient-ascent-based
		NC05 GCa10 g GCb10 g TP09 QS08 SLIC12 ERGC
		[18]	[21]	[21]	[12]	[20]	[1]
	BR	0.68	0.69	0.70	0.61	0.79	0.82	0.88
	UE	0.22	0.22	0.22	0.24	0.20	0.19	0.19
	ASA	N/A	N/A	N/A	N/A	N/A	97.94	98.22
	Segmentation speed
	320 × 240 image	×494	×14	×11	×22	×13	×1	×0.33
	2048 × 1536 image	N/A	×21	×15	×53	×12	×1	×0.49
	Control over amount of superpixels	Yes	Yes	Yes	Yes	No	Yes	Yes
	Control over superpixels compactness	No	No	No	No	No	Yes	Yes
	Supervoxel extension	No	Yes	Yes	No	No	Yes	Yes

Preprint submitted to IRBMDecember 12, 2013

Acknowledgements

This work is part of the SIRTDose project funded by the grant "Physique Cancer" 2012 (INSERM -Plan Cancer). The authors would like to thank the authors of [1] for the enlightenments on their algorithm.


Yves Achdou
email: achdou@ljll-univ-paris-diderot.fr

Alessio Porretta

CONVERGENCE OF A FINITE DIFFERENCE SCHEME TO WEAK SOLUTIONS OF THE SYSTEM OF PARTIAL DIFFERENTIAL EQUATION ARISING IN MEAN FIELD GAMES

Keywords: Mean field games, finite difference schemes, convergence AMS subject classifications. 65M06, 65M012, 9108, 91A23, 49L25

Mean field type models describing the limiting behavior of stochastic differential games as the number of players tends to +∞, have been recently introduced by J-M. Lasry and P-L. Lions. Under suitable assumptions, they lead to a system of two coupled partial differential equations, a forward Bellman equation and a backward Fokker-Planck equations. Finite difference schemes for the approximation of such systems have been proposed in previous works. Here, we prove the convergence of these schemes towards a weak solution of the system of partial differential equations.

1. Introduction. Mean field type models describing the asymptotic behavior of stochastic differential games (Nash equilibria) as the number of players tends to +∞ have recently been introduced by J-M. Lasry and P-L. Lions [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 termed mean field games by the same authors. Related ideas have been developed independently in the engineering literature by Huang-Caines-Malhamé, see for example [START_REF] Huang | Large-population cost-coupled LQG problems with nonuniform agents: individual-mass behavior and decentralized -Nash equilibria[END_REF]. For brevity, the acronym MFG will sometimes be used for mean field games. Examples of MFG models with applications in economics and social sciences are proposed in [START_REF] Guéant | Mean field games and applications[END_REF][START_REF] Achdou | Partial differential equation models in macroeconomics[END_REF].

The simplest MFG model lead to systems of evolutive partial differential equations involving two unknown scalar functions: the density of the agents in a given state x ∈ R d , namely m = m(t, x) and the value function u = u(t, x). The present work is devoted to finite difference schemes for the systems of partial differential equations. Although the methods and the theoretical results obtained below can be easily generalized, the present work focuses on the two-dimensional case for the following reasons: 1) the one dimensional case is easier and allows too special arguments; 2) in dimension two, the description of the discrete methods discussed below remain fairly simple. Besides, several important applications of the mean field games theory are two-dimensional, in particular those related to crowd dynamics.

In the state-periodic setting, typical MFG model comprises the following system of partial differential equations in (0, T ) × T 2 ∂u ∂t (t, x) -ν∆u(t, x) + H(x, ∇u(t, x)) = F (m(t, x)),

∂m ∂t (t, x) + ν∆m(t, x) + div m(t, x) ∂H ∂p (x, ∇u(t, x)) = 0, (

with the initial and terminal conditions u(0, x) = u 0 (x), m(T, x) = m T (x), in T 2 , (

given a cost function u 0 and a probability density m T .

Here, we denote by T 2 = [0, 1] 2 the 2-dimensional unit torus, and by ∆, ∇ and div, respectively, the Laplace, the gradient and the divergence operator acting on the state variable x. The parameter ν is the diffusion coefficient. Hereafter, we will always assume that ν > 0. The system also involves the scalar Hamiltonian H(x, p), which is assumed to be continuous, convex and C 1 regular with respect to p. The notation ∂H ∂p (x, q) is used for the gradient of p → H(x, p) at p = q. Finally, in the term F (m(t, x)), F is a continuous real valued function defined on R + . Hereafter the notation Q will be used for the space-time cylinder (0, T ) × T 2 .

We have chosen to focus on the case when the cost u 0 depends directly on x. In some realistic situations, the final cost may depend on the density of the players, i.e. u |t=0 = Φ 0 [m |t=0 ](x), where Φ 0 is an operator acting on probability densities, which may be local or not. We will not tackle this aspect, in order to keep the discussion as simple as possible. Similarly, by working on the torus T 2 , we avoid the discussion of the boundary conditions, but other boundary value problems with for example Dirichlet or Neumann conditions could be considered. It is also possible to consider different initial conditions than in (1.3): if there is a condition of the type m(t = 0, •) = m 0 instead of u(t = 0, •) = u 0 , then the system models a planning problem, see [START_REF] Lions | Cours du Collège de France[END_REF] for a description of the model and mathematical results, and [START_REF] Porretta | On the planning problem for a class of mean field games[END_REF][START_REF] Porretta | On the planning problem for the mean field games system[END_REF] for new existence and uniqueness results. System (1.1)-(1.2) consists then of a forward Bellman equation coupled with a backward Fokker-Planck equation. The forward-backward structure is an important feature of this system, which makes it necessary to design new strategies for its mathematical analysis (see [START_REF] Lasry | Jeux à champ moyen. II. Horizon fini et contrôle optimal[END_REF][START_REF] Lasry | Mean field games[END_REF]) and for numerical approximation. The main results on the mathematical analysis of (1.1)-(1.2) are contained in the pioneering articles [START_REF] Lasry | Jeux à champ moyen. II. Horizon fini et contrôle optimal[END_REF][START_REF] Lasry | Mean field games[END_REF], but many important aspects of the theory developed by J-M. Lasry and P-L. Lions on MFG are not published in journals or books. They can nevertheless be found in the videos of the lectures of P-L. Lions (in French) at Collège de France: see [START_REF] Lions | Cours du Collège de France[END_REF]. A very good introduction is also given in the notes by P. Cardaliaguet, [START_REF] Cardaliaguet | Notes on mean field games[END_REF], with a special emphasis on the deterministic case, i.e. ν = 0 in (1.1)-(1.2). The survey of Gomes et al [START_REF] Gomes | Mean field games models-a brief survey[END_REF] also addresses interesting extensions of the model, and the so-called master equation first introduced in [START_REF] Lions | Cours du Collège de France[END_REF].

Depending on the data and on F and H, different notions of solutions can be relevant for (1.1)-(1.3): indeed, if the right hand side of (1.1) is replaced by Φ[m(t, •)](x) where Φ is a nonlocal smoothing operator, mapping probability measures on T 2 to C 1 functions, if H depends smoothly on x and if the data u 0 and m T are smooth, then classical solutions can be found, see [START_REF] Lasry | Jeux à champ moyen. II. Horizon fini et contrôle optimal[END_REF][START_REF] Lasry | Mean field games[END_REF]. The same is true if e.g. H is Lipschitz continuous w.r.t. its second argument p and F in (1.1) is a continuous function. The situation is different in the case when H has a strictly superlinear growth with respect to p and F is a continuous function: in this case, one has to look for weak solutions, see [START_REF] Lasry | Mean field games[END_REF] and the recent article [START_REF] Porretta | Weak Solutions to Fokker-Planck Equations and Mean Field Games[END_REF] which is devoted to weak solutions to Fokker-Planck equations and to the system (1.1)-(1.2). Since the (semi-)analytic solutions of the MFG system do not exist in general, any attempt to apply MFG models and to get qualitative/quantitative information from them must rely on numerical simulations and scientific computing. Therefore, the research has also been active on numerical methods for approximating (1.1)-(1.3): a numerical method based on the reformulation of the model as an optimal control problem for the Fokker-Planck equation with an application in economics was proposed in [START_REF] Lachapelle | Computation of mean field equilibria in economics[END_REF]. Discrete time, finite state space mean field games were discussed in [START_REF] Gomes | Discrete time, finite state space mean field games[END_REF]. We also refer to [START_REF] Guéant | Mean field games equations with quadratic Hamiltonian: a specific approach[END_REF][START_REF] Guéant | New numerical methods for mean field games with quadratic costs[END_REF] for a specific constructive approach when the Hamiltonian is quadratic. Finally, semi-Lagrangian approximations have been studied in [START_REF] Carlini | A fully discrete semi-Lagrangian scheme for a first order mean field game problem[END_REF][START_REF] Carlini | A semi-lagrangian scheme for a degenerate second order mean field game system[END_REF]. The finite difference method described and studied below has first been proposed and discussed in [START_REF] Achdou | Mean field games: numerical methods[END_REF][START_REF] Achdou | Mean field games: numerical methods for the planning problem[END_REF]. It will be reviewed in § 2. The numerical scheme basically relies on monotone approximations of the Hamiltonian and on a suitable weak formulation of the Fokker-Planck equation. It has several important features:

• existence and uniqueness for the discretized problems can be obtained by similar arguments as those used in the continuous case,

• they are robust when ν → 0 (the deterministic limit of the models),

• bounds on the solutions, which are uniform in the grid step, can be proved under reasonable assumptions on the data. A first result on the convergence to classical solutions was given in [START_REF] Achdou | Mean field games: numerical methods[END_REF]. The issue of convergence was studied with more details in [START_REF] Achdou | Mean field games: convergence of a finite difference method[END_REF][START_REF] Achdou | Finite difference methods for mean field games[END_REF]: in these works, the starting point/assumption was the existence of a classical solution of (1.1)- (1.3). The proof of convergence mainly consisted in plugging the classical solution into the system of equations arising from the finite difference method, and use the consistency and stability properties of the scheme in order to get estimates and pass to the limit. In the present work, the goal is different: we wish to prove that as the grid steps tend to zero, the solution of the discretized MFG system converges to a weak solution of (1.1)-(1.3), without assuming the existence of the latter; so this work will supply as a by-product a new strategy for proving the existence of weak solutions. One key step will be to obtain a priori estimates on the solutions to the discrete systems, and these will mainly come from the fact that the structure of the MFG system is preserved by the chosen finite difference method. This step will be achieved in § 4. Note that § 4.3 is concerned with a priori estimates for the discrete version of the Fokker-Planck equation (1.2). These estimates may have their own interest, independently from MFG models, and may be put in relation with recent works of Gallouët et al [START_REF] Gallouët | Convergence of a finite volume scheme for the convection-diffusion equation with L 1 data[END_REF] in the context of finite volume methods. Once these estimates are obtained, the most important difficulty will be to pass to the limit in the discrete Bellman equation. The strategy to that purpose is to first prove some L 1 compactness of the sequence of state-gradients of the discrete solutions, then to adapt some techniques proposed by Boccardo, Murat and Puel, see [START_REF] Boccardo | Existence de solutions non bornées pour certaines équations quasi-linéaires[END_REF], for studying weak solutions of (1.1). This will done in § 5.1 and 6. Here also, we think that the passage to the limit in the discrete Bellman equation may have an interest for itself.

2. Finite difference schemes. In the present paragraph, we discuss the finite difference method originally proposed in [START_REF] Achdou | Mean field games: numerical methods[END_REF]. Let N T be a positive integer and ∆t = T /N T , t n = n∆t, n = 0, . . . , N T . Let T 2 h be a uniform grid on the torus with mesh step h, (assuming that 1/h is an integer N h ), and x ij denote a generic point in T 2

h . The values of u and m at (x i,j , t n ) are respectively approximated by u n i,j and m n i,j . Let u n (resp. m n ) be the vector containing the values u n i,j (resp. m n i,j ), for 0 ≤ i, j < N h indexed in the lexicographic order. Hereafter, such vectors will be termed grid functions on T 2 h or simply grid functions. For all grid functions z, all i and j, we agree that z i,j = z (i mod N h ),(j mod N h ) .

Elementary finite difference operators. Let us introduce the elementary finite difference operators

(D + 1 u) i,j = u i+1,j -u i,j h and (D + 2 u) i,j = u i,j+1 -u i,j h , (2.1)
and define D h u as the grid function with values in R 2 :

(D h u) i,j = (D + 1 u) i,j , (D + 2 u) i,j ∈ R 2 . (2.2)
Let [∇ h u] i,j be the collection of the four possible one sided finite differences at x i,j :

[∇ h u] i,j = (D + 1 u) i,j , (D + 1 u) i-1,j , (D + 2 u) i,j , (D + 2 u) i,j-1 ∈ R 4 . (2.3)
We will also need the standard five point discrete Laplace operator

(∆ h u) i,j = - 1 h 2 (4u i,j -u i+1,j -u i-1,j -u i,j+1 -u i,j-1 ).
For a set v = (v n ) n=0,...,N T , where v n is grid functions on T 2 h , it will be convenient to define the family of grid functions:

∂ t,∆t v ≡ v n+1 -v n ∆t n=0,...,N T -1
.

(2.4)

Numerical Hamiltonian. In order to approximate the term H(x, ∇u) in (1.1), we consider a numerical Hamiltonian g : T 2 ×R 4 → R, (x, q 1 , q 2 , q 3 , q 4 ) → g (x, q 1 , q 2 , q 3 , q 4 ). Hereafter we will often assume that the following conditions hold:

(g 1 ) monotonicity: g is nonincreasing with respect to q 1 and q 3 and nondecreasing with respect to q 2 and q 4 . (g 2 ) consistency: g (x, q 1 , q 1 , q 2 , q 2 ) = H(x, q), ∀x ∈ T 2 , ∀q = (q 1 , q 2 ) ∈ R 2 . (g 3 ) regularity: g is continuous and of class C 1 w.r.t. (q 1 , q 2 , q 3 , q 4 ). (g 4 ) convexity : (q 1 , q 2 , q 3 , q 4 ) → g (x, q 1 , q 2 , q 3 , q 4 ) is convex. We will approximate H(•, ∇u)(x i,j ) by g(x i,j , [∇ h u] i,j ). Standard examples of numerical Hamiltonians fulfilling these requirements are provided by Lax-Friedrichs or upwind schemes, see [START_REF] Achdou | Mean field games: numerical methods[END_REF]. For Hamiltonians of the form

H(x, p) = H(x) + |p| β , β ∈ (1, ∞), we may choose g(x, q) = H(x) + G(q - 1 , q + 2 , q - 3 , q + 4 ), (2.5)
where, for a real number r, r + = max(r, 0) and r -= max(-r, 0) and where G : (R + ) 4 → R + is given by

G(p) = |p| β = (p 2 1 + p 2 2 + p 2 3 + p 2 4 ) β 2 .
(2.6)

Discrete Bellman equation. The discrete version of the Bellman equation is obtained by applying a semi-implicit Euler scheme to (1.1),

u n+1 i,j -u n i,j ∆t -ν(∆ h u n+1 ) i,j + g(x i,j , ∇ h u n+1 i,j ) = F (m n i,j ), (2.7)
for all points in T 2 h and all n, 0 ≤ n < N T , where all the discrete operators have been introduced above. Given (m n ) n=0,...,N T -1 , (2.7) and the initial condition u 0 i,j = u 0 (x i,j ) for all (i, j) completely characterizes (u n ) 0≤n≤N T .

Discrete Fokker-Planck equation. In order to approximate equation (1.2), it is convenient to consider its weak formulation which involves in particular the term

T 2 div m ∂H ∂p (x, ∇u) w(x) dx.
By periodicity,

T 2 div m ∂H ∂p (x, ∇u) w(x) dx = - T 2 m(x) ∂H ∂p (x, ∇u(x)) • ∇w(x) dx
holds for any test function w. The right hand side in the identity above will be approximated by

-h 2 i,j m i,j ∇ q g(x i,j , [∇ h u] i,j ) • [∇ h w] i,j = h 2 i,j T i,j (u, m)w i,j ,
where the transport operator T is defined as follows:

T i,j (u, m) = 1 h                 m i,j ∂g ∂q 1 (x i,j , [∇ h u] i,j ) -m i-1,j ∂g ∂q 1 (x i-1,j , [∇ h u] i-1,j ) +m i+1,j ∂g ∂q 2 (x i+1,j , [∇ h u] i+1,j ) -m i,j ∂g ∂q 2 (x i,j , [∇ h u] i,j )    +    m i,j ∂g ∂q 3 (x i,j , [∇ h u] i,j ) -m i,j-1 ∂g ∂q 3 (x i,j-1 , [∇ h u] i,j-1 ) +m i,j+1 ∂g ∂q 4 (x i,j+1 , [∇ h u] i,j+1 ) -m i,j ∂g ∂q 4 (x i,j , [∇ h u] i,j )                 . (2.8)
The discrete version of equation (1.2) is chosen as follows:

m n+1 i,j -m n i,j ∆t + ν(∆ h m n ) i,j + T i,j (u n+1 , m n ) = 0, (2.9)
for all n = 0, . . . , N T -1. This scheme is implicit w.r.t. to m and explicit w.r.t. u because the considered Fokker-Planck equation is backward. Given u this is a system of linear equations for m. We introduce the compact and convex set

K h = {(m i,j ) 0≤i,j<N h : h 2 i,j m i,j = 1; m i,j ≥ 0} (2.10)
which can be viewed as the set of the discrete probability measures. It is easy to see that if m n satisfies (2.9) for 0

≤ n < N T and if m N T ∈ K h , then m n ∈ K h for all n, 0 ≤ n < N T . Remark 1.
An important property of T is that the operator m → -ν(∆ h m) i,j -T i,j (u, m) i,j is the adjoint of the linearized version of the operator u → -ν(∆ h u) i,j + g(x i,j , [∇ h u] i,j ) i,j . This property implies that the structure of (1.1)-(1.2) is preserved in the discrete version (2.7)-(2.9). In particular, it implies the uniqueness result stated in Theorem 2.2 below.

Summary. The fully discrete scheme for system (1.1),(1.2),(1.3) is therefore the following: for all 0 ≤ i, j < N h and 0

≤ k < N T    u k+1 i,j -u k i,j ∆t -ν(∆ h u k+1 ) i,j + g(x i,j , ∇ h u k+1 i,j ) = F (m k i,j ), m k+1 i,j -m k i,j ∆t + ν(∆ h m k ) i,j + T i,j (u k+1 , m k ) = 0, (2.11)
with the initial and terminal conditions

u 0 i,j = u 0 (x i,j ), m N T i,j = 1 h 2 |x-xi,j |∞≤h/2 m T (x)dx, 0 ≤ i, j < N h . (2.
12)

The following theorem was proved in [START_REF] Achdou | Mean field games: numerical methods[END_REF] (using essentially Brouwer's fixed point theorem and estimates on the solutions of the discrete Bellman equation): Theorem 2.1. Assume that (g 1 )-(g 3 ) hold, that u 0 is continuous on T 2 and that m T ∈ L 1 (T 2 ) is a probability density, i.e. m T ≥ 0 and T 2 m T (x)dx = 1; then (2.11)-(2.12) has a solution such that m n ∈ K h , ∀n.

Since (2.11)-(2.12) has exactly the same structure as the continuous problem (1.1)-(1.3), uniqueness has been obtained in [START_REF] Achdou | Mean field games: numerical methods[END_REF] with the same arguments as in [START_REF] Lasry | Jeux à champ moyen. II. Horizon fini et contrôle optimal[END_REF]:

Theorem 2.2. Assume that (g 1 )-(g 4 ) hold and that F is nondecreasing then (2.11)-(2.12) has a unique solution.

Remark 2. Efficient algorithms for solving system (2.11)-(2.12) require special efforts, essentially because of the forward-backward structure already discussed above. We refer to [START_REF] Achdou | Mean field games: numerical methods[END_REF] for the description of possible algorithms and numerical results.

3. Running assumptions and statement of the main result. We now summarize the assumptions that will be made in the whole work.

• u 0 is a continuous function on T 2 • m T is a nonnegative function in L ∞ (T 2 ) such that T 2 m T (x)dx = 1 • F is a continuous function on R + , which is bounded from below.
• The Hamiltonian1 (x, p) → H(x, p) is assumed to be convex with respect to p and C 1 regular w.r.t. x and p. • The discrete Hamiltonian g satisfies (g 1 )-(g 4 ) and the further assumption (g 5 ) There exist positive constants c 1 , c 2 , c 3 , c 4 such that

g q (x, q) • q -g(x, q) ≥ c 1 |g q (x, q)| 2 -c 2 , (3.1
)

|g q (x, q)| ≤ c 3 |q| + c 4 . (3.2)
Take for example g as in (2.5) (2.6). It is clear that

g q (x, q) • q = βG(q - 1 , q + 2 , q - 3 , q + 4 ), hence g q (x, q) • q -g(x, q) = (β -1)G(q - 1 , q + 2 , q - 3 , q + 4 ) -H(x). Since |g q (x, q)| 2 = β 2 G(q - 1 , q + 2 , q - 3 , q + 4 ) 2 β-1 β , we see that (g 5 ) holds if 1 ≤ β ≤ 2.
We can now state the main result of this article, which establishes the convergence of the solutions of the finite difference scheme towards a weak solution of the continuous mean field games system. Theorem 3.1. Let (u n ), (m n ) be a solution of the discrete system (2.11)-(2.12) and u h,∆t , m h,∆t be the piecewise constant functions which take the values u n+1 i,j and

m n i,j , respectively, in (t n , t n+1 ) × (ih -h/2, ih + h/2) × (jh -h/2, jh + h/2).
There exists a subsequence of h and ∆t (not relabeled) and functions ũ, m, which belong to L α (0, T ; W 1,α (T 2 )) for any α ∈ [1, 4 3 ), such that u h,∆t → ũ and m h,∆t → m in L β (Q) for all β ∈ [1, 2), and (ũ, m) is a weak solution to the system (1.1)-(1.3) in the following sense:

(i) H(•, Dũ) ∈ L 1 (Q), mF ( m) ∈ L 1 (Q), m[H p (•, Dũ)•Dũ-H(•, Dũ)] ∈ L 1 (Q) (ii) (ũ, m) satisfies (1.1)-(1.2) in the sense of distributions (iii) ũ, m ∈ C 0 ([0, T ]; L 1 (T 2 )) and ũ| t=0 = u 0 , m| t=T = m T .
Remark 3. We recall that, if F is nondecreasing and p → H(x, p) is strictly convex at infinity, it is proved in [START_REF] Porretta | Weak Solutions to Fokker-Planck Equations and Mean Field Games[END_REF] that weak solutions are unique whenever H satisfies the structure conditions

H p (t, x, p) • p ≥ r H(t, x, p) -γ |H p (t, x, p)| ≤ β (1 + |p| r-1 ) H(t, x, p) ≥ α|p| r -γ
for some r ∈ (1, 2] and some positive constant α, β, γ.

Therefore, in this case the convergence established in the above theorem holds for the whole sequence, and not only for a subsequence.

A priori estimates.

4.1. Norms and semi-norms. It is useful to define the following norms and semi-norms: for a grid function v ≡ (v i,j ) i,j , we define

v L s (T 2 h ) =   h 2 i,j |v i,j | s   1 s , (4.1
)

|v| W 1,s (T 2 h ) =   h 2 i,j (D + 1 v i,j ) 2 + (D + 2 v i,j ) 2 s 2   1 s (4.2) v W 1,s (T 2 h ) = v s L s (T 2 h ) + |v| s W 1,s (T 2 h ) 1 s , (4.3)
where D + 1 v and D + 2 v are defined in (2.1). We shall also write |v|

H 1 (T 2 h ) = |v| W 1,2 (T 2 h ) , v H 1 (T 2 h ) = v W 1,2 (T 2 h )
, and define the discrete L 2 scalar product:

(v, w) L 2 (T 2 h ) = h 2 i,j v i,j w i,j .
We recall the discrete Sobolev inequality: for any s < ∞, there exists a constant C such that for any grid function v,

v L s (T 2 h ) ≤ C v L 2 (T 2 h ) + |v| H 1 (T 2 h ) . For s > 1, we define the dual norm v W -1,s (T 2 h ) , 1 s + 1 s = 1 by v W -1,s (T 2 h ) = sup w =0 (v, w) L 2 (T 2 h ) w W 1,s (T 2 h )
.

Define Q h,∆t = ∆t{0, . . . , N T -1} × T 2 h .
For a function w defined on Q h,∆t , w ≡ (w n i,j ) i,j , 0 ≤ n ≤ N T , we define for s ∈ [1, +∞),

w L s (Q h,∆t ) = ∆t N T n=0 w n s L s (T 2 h ) 1 s . (4.4)
4.2. First estimates. Hereafter, the constants appearing in the a priori estimates, for example c, C, are independent of h and ∆t. In this paragraph, we state the first a priori estimates stemming from the structure of the system. Although we have already given the set of running assumptions, we think that it may be useful to specify which assumptions are really required by each particular result.

Lemma 4.1. Under Assumptions (g 1 ) and (g 3 ), if F is bounded from below by a constant F , u 0 is continuous on T 2 , then for all i, j, n,

u n i,j ≥ u -T F -max x∈T 2 H(x, 0) - , (4.5)
where u = min x∈T 2 u 0 (x). Lemma 4.2. Under Assumptions (g 1 ), (g 3 ) and (g 5 ), if F is bounded from below by F , u 0 is continuous on T 2 and m T is bounded from above by mT , then there exists a constant C such that

h 2 ∆t N T -1 k=0 i,j m k i,j g q (x i,j , ∇ h u k+1 i,j ) 2 ≤ C, (4.6)
h 2 ∆t N T -1 k=0 i,j g(x i,j , ∇ h u k+1 i,j ) ≤ C, (4.7)
h 2 ∆t N T -1 k=0 i,j m k i,j F (m k i,j ) ≤ C. (4.8)
Proof. Consider ũn i,j = n∆tF ( mT ) for all i, j, n. We get immediately

ũn+1 i,j -ũn i,j ∆t -ν(∆ h ũn+1 ) i,j = F ( mT ) (4.9)
Subtract (4.9) from (2.7) and multiply the resulting equation by m n i,j -mT . Similarly, multiply (2.9) by u n+1 i,j -ũn+1 i,j . Adding the two resulting identities and summing with respect to n, one gets:

h 2 ∆t N T -1 k=0 i,j m k i,j g q (x i,j , ∇ h u k+1 i,j ) • ∇ h u k+1 i,j -g(x i,j , ∇ h u k+1 i,j ) + h 2 ∆t N T -1 k=0 i,j mT g(x i,j , ∇ h u k+1 i,j ) + h 2 ∆t N T -1 k=0 i,j (m k i,j -mT )(F (m k i,j ) -F ( mT )) =(m N T -mT , u N T -T F ( mT )) L 2 (T 2 h ) -(m 0 -mT , u 0 ) L 2 (T 2 h ) . (4.10)
1. Since m N T -mT is nonpositive with a bounded mass, and since u n is bounded from below, see (4.5), the term (m

N T -mT , u N T -T F ( mT )) L 2 (T 2 h )
in the right hand side of (4.10) is bounded from above by a constant independent of h and ∆t. 2. It is straightforward to see that (m 0 -mT , u 0 ) L 2 (T 2 h ) ≤ (1 + mT ) u 0 ∞ . 3. Since F is continuous, there exists a constant c such that F (t) ≤ 1 2 mT tF (t)+ c, ∀t ≥ 0. Hence,

h 2 ∆t N T -1 k=0 i,j (m k i,j -mT )F (m k i,j ) ≥ 1 2 h 2 ∆t N T -1 k=0 i,j m k i,j F (m k i,j ) -c. 4. Finally, h 2 ∆t N T -1 k=0 i,j (m k i,j -mT )F ( mT ) = T (1 -mT )F (

A priori estimates from the discrete Fokker-Planck equation.

The following estimates for the Fokker-Planck equation may have their own interest: Lemma 4.3. Assume (g 1 ) and (g 3 ). Let ψ be a non decreasing and concave function defined on R + . For any grid function v = (v i,j ), any positive grid function m = (m i,j ) and any positive number η, i,j

T i,j (v, m)ψ(m i,j ) ≤ η 2 i,j [∇ h ψ(m)] i,j • [∇ h m] i,j + 1 2η i,j m 2 i,j ψ (m i,j ) ∂g ∂q 1 (x i,j , [∇ h v] i,j ) 2 1 {mi+1,j >mi,j } + 1 2η i,j m 2 i,j ψ (m i,j ) ∂g ∂q 2 (x i,j , [∇ h v] i,j ) 2 1 {mi-1,j >mi,j } + 1 2η i,j m 2 i,j ψ (m i,j ) ∂g ∂q 3 (x i,j , [∇ h v] i,j ) 2 1 {mi,j+1>mi,j } + 1 2η i,j m 2 i,j ψ (m i,j ) ∂g ∂q 4 (x i,j , [∇ h v] i,j ) 2 1 {mi,j-1>mi,j } . (4.11)
In particular, if m does not vanish, then for ψ(z) = ln(z), i,j

T i,j (v, m) ln(m i,j ) ≤ η 2 i,j [∇ h ln(m)] i,j • [∇ h m] i,j + 1 2η i,j m i,j |g q (x i,j , [∇ h v] i,j )| 2 . (4.12)
Proof. By the definition of T , we can split the sum S = i,j T i,j (v, m)ψ(m i,j ) as follows:

S = - i,j m i,j ∇ q g(x i,j , [∇ h v] i,j ) • [∇ h ψ(m)] i,j = S 1 + S 2 + S 3 + S 4 ,
where

S 1 = - 1 h i,j m i,j ∂g ∂q 1 (x i,j , [∇ h v] i,j )(ψ(m i+1,j ) -ψ(m i,j )), S 2 = - 1 h i,j m i,j ∂g ∂q 2 (x i,j , [∇ h v] i,j )(ψ(m i,j ) -ψ(m i-1,j )), S 3 = - 1 h i,j m i,j ∂g ∂q 3 (x i,j , [∇ h v] i,j )(ψ(m i,j+1 ) -ψ(m i,j )), S 4 = - 1 h i,j m i,j ∂g ∂q 4 (x i,j , [∇ h v] i,j )(ψ(m i,j ) -ψ(m i,j-1 )).
It is enough to focus on S 1 since the same arguments can be used for the other sums. Since g is nonincreasing w.r.t. q 1 ,

S 1 ≤ - 1 h i,j m i,j ∂g ∂q 1 (x i,j , [∇ h v] i,j )(ψ(m i+1,j ) -ψ(m i,j )) + .
Since ψ is nondecreasing, if m i+1,j > m i,j , the factor (ψ(m i+1,j ) -ψ(m i,j )) + can be rewritten

ψ(m i+1,j ) -ψ(m i,j ) m i+1,j -m i,j 1 2 (ψ(m i+1,j ) -ψ(m i,j ))(m i+1,j -m i,j ) 1 2 .
Since ψ is nondecreasing and concave, m i+1,j > m i,j implies that

0 ≤ ψ(m i+1,j ) -ψ(m i,j ) m i+1,j -m i,j ≤ ψ (m i,j ).
Hence, if m i+1,j ≥ m i,j , then

(ψ(m i+1,j ) -ψ(m i,j )) + ≤ ψ (m i,j )(ψ(m i+1,j ) -ψ(m i,j ))(m i+1,j -m i,j ) 1 2
, which implies that

S 1 ≤ 1 2η i,j m 2 i,j ψ (m i,j ) ∂g ∂q 1 (x i,j , [∇ h v] i,j ) 2 1 {mi+1,j >mi,j } + η 2h 2 i,j
(ψ(m i+1,j ) -ψ(m i,j ))(m i+1,j -m i,j )1 {mi+1,j >mi,j } .

Lemma 4.4. Assume (g 1 ) and (g 3 ). Let ψ be a non decreasing and concave function defined on R + . For any positive grid functions (m k i,j ), k = 0, . . . , N T ,

N T k=n+1 i,j m k i,j (ψ(m k i,j ) -ψ(m k-1 i,j )) ≥ - N T k=n+1 i,j m k i,j ψ (m k i,j )(m k-1 i,j -m k i,j ) + .
(4.13) If m k ∈ K h for all k ∈ {0, . . . , N T } and does not vanish, then

h 2 N T -1 k=n i,j m k+1 i,j -m k i,j ln(m k i,j ) ≤ h 2 i,j m N T i,j ln(m N T i,j ) -h 2 i,j
m n i,j ln(m n i,j ) + 1.

(4.14)

Proof. Since ψ is non decreasing,

N T k=n+1 i,j m k i,j (ψ(m k i,j )-ψ(m k-1 i,j )) ≥ N T k=n+1 i,j m k i,j (ψ(m k i,j )-ψ(m k-1 i,j ))1 {m k i,j <m k-1 i,j } .
From the concavity of ψ, if

m k i,j < m k-1 i,j , then ψ(m k i,j )-ψ(m k-1 i,j ) ≥ ψ (m k i,j )(m k i,j - m k-1 i,j
). Then (4.13) follows from the last two points. Let us turn to (4.14): for any > 0,

h 2 N T -1 k=n i,j m k+1 i,j -m k i,j ln(m k i,j + ) =h 2 i,j m N T i,j ln(m N T i,j + ) -h 2 i,j m n i,j ln(m n i,j + ) -h 2 N T k=n+1 i,j m k i,j (ln(m k i,j + ) -ln(m k-1 i,j + )),
and (4.13) with ψ(z) = ln(z + ) yields

-h 2 N T k=n+1 i,j m k i,j (ln(m k i,j + ) -ln(m k-1 i,j + )) ≤ h 2 N T k=n+1 i,j m k i,j m k i,j + (m k-1 i,j -m k i,j ) + ≤ h 2 N T k=n+1 i,j (m k-1 i,j -m k i,j ) + ≤ 1,
where the last estimate comes from the fact that the grid functions m k all belong to K h . Hence,

h 2 N T -1 k=n i,j m k+1 i,j -m k i,j ln(m k i,j + ) ≤ h 2 i,j m N T i,j ln(m N T i,j + )-h 2 i,j
m n i,j ln(m n i,j + )+1. and (4.14) is obtained by letting tend to 0. Lemma 4.5. If m N T ∈ K h and (g 1 ) (g 3 ) hold, then there exists a constant C such that, for any number η, 0 < η < ν, a solution (m n i,j ) of (2.9) satisfies

max n h 2 i,j m n i,j | ln(m n i,j )| + (ν -η)∆t N T -1 k=0 √ m k 2 H 1 (T 2 h ) ≤ C + h 2 i,j m N T i,j | ln(m N T i,j )| + h 2 ∆t 2η N T -1 k=0 i,j m k i,j g q x i,j , ∇ h u k+1 i,j 2 . (4.15)
For all α ∈ [1, 2), there exists a constant c such that

m α L α (Q h,∆t ) ≤c   1 + h 2 i,j m N T i,j | ln(m N T i,j )| + h 2 ∆t N T -1 k=0 i,j m k i,j g q (x i,j , ∇ h u k+1 i,j ) 2   . (4.16)
Proof.

Step 1. Take > 0 and consider mn i,j = m n i,j + . Note that mn i,j > 0 for all i, j, n. Multiply the second equation of (2.11) by ln( mn i,j ) and sum for all i, j and k = n, . . . , N T -1:

0 =h 2 N T -1 k=n i,j mk+1 i,j -mk i,j ln( mk i,j ) - νh 2 ∆t 2 N T -1 k=n i,j ∇ h mk i,j • ∇ h ln( mk ) i,j + h 2 ∆t N T -1 k=n i,j T i,j (u k+1 , m k ) ln( mk i,j ).
From (4.12) and (4.14), we deduce that

h 2 i,j mn i,j ln( mn i,j ) + (ν -η) h 2 ∆t 2 N T -1 k=n i,j ∇ h mk i,j • ∇ h ln( mk ) i,j ≤1 + h 2 i,j mN T i,j ln( mN T i,j ) + h 2 ∆t 2η N T -1 k=n i,j mk i,j g q (x i,j , ∇ h u k+1 i,j ) 2 ,
and since mk i,j ln( mk i,j ) ≥ -e -1 , max

n h 2 i,j mn i,j | ln( mn i,j )| + ν -η 2 h 2 ∆t N T -1 k=0 i,j ∇ h mk i,j • ∇ h ln( mk ) i,j ≤ C + h 2 i,j mN T i,j | ln( mN T i,j )| + h 2 ∆t 2η N T -1 k=0 i,j mk i,j g q (x i,j , ∇ h u k+1 i,j ) 2 .
(4.17) Consider now the quantity

| √ mk | 2 H 1 (T 2 h ) , i.e. | √ mk | 2 H 1 (T 2 h ) = i,j mk i+1,j -mk i,j 2 + i,j mk i,j+1 -mk i,j 2 .
Since mk i,j > 0, we can write mk i+1,j -mk

i,j 2 = mk i,j 1 + h (D + 1 mk )i,j mk i,j - 1
2
where (D + 1 mk ) i,j is defined in (2.1). Since the inequality ( √ 1 + z -1) 2 ≤ z ln(1 + z) holds for any number z ≥ -1, we infer that

mk i+1,j -mk i,j 2 ≤ h(D + 1 mk ) i,j ln(1 + h (D + 1 m) i,j mk i,j ) = h(D + 1 mk ) i,j ln( mk i+1,j ) -ln( mk i,j ) = h 2 (D + 1 mk ) i,j (D + 1 ln( mk )) i,j .
Since the same kind of estimate holds for mk i,j+1 -mk i,j

2

, we obtain that

| √ mk | 2 H 1 (T 2 h ) ≤ h 2 i,j D h mk i,j • D h ln( mk ) i,j , (4.18)
and the fact that m satisfies (4.15) follows from (4.17) and (4.18). Let us now prove (4.16): consider α ∈ [1, 2): there exists a unique number p ≥ 1 such that 1 α = 1 2 + 1 2p : we have the interpolation inequality mk

L α (T 2 h ) ≤ mk 1 2 L 1 (T 2 h ) mk 1 2 L p (T 2 h ) . But mk L 1 (T 2 h ) = 1 + and mk 1 2 L p (T 2 h ) = √ mk L 2p (T 2
h ) . From the discrete Sobolev inequalities, we deduce that

mk α L α (T 2 h ) ≤ (1 + ) α 2 √ mk α L 2p (T 2 h ) ≤ C √ mk α H 1 (T 2 h ) ≤ C 1 + √ mk 2 H 1 (T 2 h ) ,
which yields that m satisfies (4.16) by summing for all k and using (4.15).

Step 2. We obtain that m satisfies (4.15) and (4.16) by letting tend to 0. Corollary 4.6. With the same assumptions as in Lemma 4.5, for any α ∈ [1, 4/3), there exists a constant c such that

D h m α L α (Q h,∆t ) + ∆t N T -1 k=0 m k+1 -m k ∆t α W -1,α (T 2 h )) ≤c   1 + h 2 i,j m N T i,j | ln(m N T i,j )| + h 2 ∆t N T -1 k=0 i,j m k i,j g q (x i,j , ∇ h u k+1 i,j ) 2   . (4.19)
Proof. Take α ∈ [1, 4/3). We start by observing that

D h m α L α (Q h,∆t ) ≤ Ch 2 ∆t N T -1 k=0 i,j |D + 1 m k | α i,j + |D + 2 m k | α i,j .
Let us estimate

N T -1 k=0 i,j |D + 1 m k | α i,j
by some quantity depending on

N T -1 k=0
i,j (D + 1 m k ) i,j (D + 1 ln(m k )) i,j : a Hölder inequality yields that

N T -1 k=0 i,j |D + 1 m k | α i,j ≤   N T -1 k=0 i,j |D + 1 m k | 2 i,j (D + 1 ln(m k )) i,j (D + 1 m k ) i,j   α 2   N T -1 k=0 i,j (D + 1 m k ) i,j (D + 1 ln(m k )) i,j α 2-α   1-α 2 =   N T -1 k=0 i,j (D + 1 m k ) i,j (D + 1 ln(m k )) i,j   α 2   N T -1 k=0 i,j (D + 1 m k ) i,j (D + 1 ln(m k )) i,j α 2-α   1-α 2 .
Standard calculus yields that

(D + 1 m k )i,j (D + 1 ln(m k ))i,j ≤ max(m k i,j , m k i+1,j ) ≤ m k i,j + m k i+1,j , therefore h 2 ∆t N T -1 k=0 i,j |D + 1 m k | α i,j ≤h 2 ∆t   N T -1 k=0 i,j (D + 1 m k ) i,j (D + 1 ln(m k )) i,j   α 2   N T -1 k=0 i,j (m k i,j + m k i+1,j ) α 2-α   1-α 2 ≤C m α 2 L α 2-α (Q h,∆t )   h 2 ∆t N T -1 k=0 i,j (D + 1 m k ) i,j (D + 1 ln(m k )) i,j   α 2
.

Note that 1 ≤ α 2-α < 2. The estimate on D h m α L α (Q h,∆t ) follows by using the same argument for bounding

N T -1 k=0 i,j |D + 2 m k | α i,j
and then (4.17) and (4.16). From (2.9), we deduce that for all grid functions w,

h 2 i,j w i,j m n+1 i,j -m n i,j ∆t =νh 2 i,j (∇ h m n ) i,j • (∇ h w) i,j + h 2 i,j m n i,j g q (x i,j , [∇ h u n+1 ] i,j ) • [∇ h w] i,j .
Hence

m n+1 -m n ∆t , w L 2 (T 2 h ) ≤ ν|m n | W 1,α (T 2 h ) |w| W 1,α (T 2 h ) +   h 2 i,j m n i,j |g q (x i,j , [∇ h u n+1 ] i,j )| 2   1 2 √ m n L 2α 2-α (T 2 h ) |w| W 1,α (T 2 h ) ,
and ∆t

N T -1 k=0 m k+1 -m k ∆t α W -1,α (T 2 h ) 1 α ≤ ν D h m L α (Q h,∆t ) +   h 2 ∆t N T -1 k=0 i,j m k i,j g q x i,j , ∇ h u k+1 i,j 2   1 2 m 1 2 L α 2-α (Q h,∆t ) . (4.20)
Note that α 2-α < 2. From (4.16),

m 1 2 L α 2-α (Q h,∆t ) ≤c   1 + h 2 i,j m N T i,j | ln(m N T i,j )| + h 2 ∆t N T -1 k=0 i,j m k i,j g q (x i,j , ∇ h u k+1 i,j ) 2   2-α 2α . (4.21)
The desired estimate on ∆t

N T -1 k=0 m k+1 -m k ∆t α W -1,α (T 2 h )
follows from (4.20)-(4.21)

and from the estimate on D h m α L α (Q h,∆t ) . Collecting the above results together with Lemma 4.2, we obtain the following conclusion:

Theorem 4.7. If F is continuous and bounded from below by a constant F , if (g 1 ), (g 3 ), (g 5 ) hold, if u 0 is continuous, then there exists a constant C such that a solution (u, m) of (2.11)-(2.12) satisfies (4.6)-(4.8), and for all α ∈ (1, 4/3),

D h m α L α (Q h,∆t ) + ∆t N T -1 k=0 m k+1 -m k ∆t α W -1,α (T 2 h )) 1 α ≤ C. (4.22)
5. L 1 -compactness results. In this section we prove the L 1 -compactness of D h u whenever the discrete heat equation has bounded L 1 data. More precisely, we assume that u = (u n ) n=0,...,N T satisfies

u n+1 i,j -u n i,j ∆t -ν(∆ h u n+1 ) i,j = f n i,j (5.1)
for all 0 ≤ i, j < N h and all n, 0 ≤ n < N T , where the data f = (f n i,j ) and the initial conditions u 0 = (u 0 i,j ) are supposed to satisfy

S u 0 L 1 (T 2 h ) + f L 1 (Q h,∆t ) ≤ c (5.2)
for some c independent of h and ∆t. In what follows, we reconstruct functions on Q from the grid functions u, and we prove the convergence of these functions as h and ∆t tend to 0, at least for subsequences. Lemma 5.1 below is concerned with piecewise constant functions built using u. It is similar to results that can be found in Gallouët et al, see e.g. [START_REF] Eymard | Finite volumes and nonlinear diffusion equations[END_REF][START_REF] Eymard | Finite volume methods[END_REF][START_REF] Gallouët | Convergence of a finite volume scheme for the convection-diffusion equation with L 1 data[END_REF] in the context of finite volume methods. Lemma 5.2 deals with approximations of the gradient with respect to x. It seems new to the best of our knowledge and may have an independent interest. Lemma 5.1. Let u h,∆t be the piecewise constant function which takes the value

u n+1 i,j in (t n , t n+1 ) × (ih -h/2, ih + h/2) × (jh -h/2, jh + h/2).
There exists a subsequence of h and ∆t (not relabeled) and a function ũ such that u h,∆t → ũ in in

L β (Q) for all β ∈ [1, 2). Moreover, ũ ∈ L α (0, T ; W 1,α (T 2 )) for any α ∈ [1, 4
3 ), and there exist a bounded Radon measure μ in Q and a bounded Radon measure μ0 in T 2 such that ũ is the unique solution of

∂ t ũ -ν∆ũ = μ in Q, ũ(0, •) = μ0 in T 2 .
(5.3)

Proof. Using the L 1 bounds on the data, we may show with the same argument as in § 4.3 that

u L β (Q h,∆t ) + ∆t N T n=1 |D h u n | α L α (T 2 h ) 1 α ≤ c
for any β ∈ [1, 2) and α ∈ [1, 4 3 ). From this estimate and (5.1), we deduce that ∆t

N T n=1 u n+1 -u n ∆t W -1,α (T 2
h ) is uniformly bounded. Recall that u h,∆t is the piecewise constant function which takes the value u n+1 i,j in (t n , t n+1 ) × (ih -h/2, ih + h/2) × (jh -h/2, jh + h/2). We can apply the discrete Aubin-Simon lemma in [START_REF] Gallouët | Convergence of a finite volume scheme for the convection-diffusion equation with L 1 data[END_REF] (Theorem 3.1): up to the extraction of a subsequence, u h,∆t converges to a function ũ in L 1 (Q), and in fact in L β (Q) for all β ∈ [1, 2). Moreover, ũ ∈ L α (0, T ; W 1,α (T 2 )) for any α ∈ [1, 4 3 ). Let f h,∆t be the piecewise constant function on which takes the values

f n i,j in (t n , t n+1 ) × (ih -h/2, ih + h/2) × (jh -h/2, jh + h/2).
Up to the extraction of a subsequence, f h∆t converges in the weak- * topology to some bounded Radon measure μ on Q. Call u 0 h the piecewise constant function on T 2 which takes the values

u 0 i,j in (ih -h/2, ih + h/2) × (jh -h/2, jh + h/2
). We may assume that u 0 h converges to a bounded measure μ0 on T 2 . In particular, testing (5.1) with smooth functions and passing to the limit, this implies that ũ satisfies

T 0 T 2 ũ (-ϕ t -ν∆ϕ) dxdt = T 0 T 2 ϕ dμ + T 2 ϕ(0)dμ 0 ,
for every ϕ ∈ C 2 (Q) such that ϕ(T ) = 0. Notice that ũ is the unique solution of the above weak formulation. We now define an approximation of Dũ from the grid function u. For a real number z, let floor(z) be the largest integer that does not exceed z, ceil(z) be the smallest integer that is not less than z. Let Du h,∆t be the piecewise constant function from Q to R 2 which takes the value

D + 1 u n+1 floor( i 2 ),ceil( j 2 ) , D + 2 u n+1 ceil( i 2 ),floor( j 2 ) in (t n , t n+1 ) × i h 2 , (i + 1) h 2 × j h 2 , (j + 1) h 2 .
More explicitly, Du h,∆t takes the value

D + 1 u n+1 i,j , D + 2 u n+1 i,j in (t n , t n+1 ) × ih, ih + h 2 × jh, jh + h 2 D + 1 u n+1 i,j , D + 2 u n+1 i,j-1 in (t n , t n+1 ) × ih, ih + h 2 × jh -h 2 , jh D + 1 u n+1 i-1,j , D + 2 u n+1 i,j-1 in (t n , t n+1 ) × ih -h 2 , ih × jh -h 2 , jh D + 1 u n+1 i-1,j , D + 2 u n+1 i,j in (t n , t n+1 ) × ih -h 2 , ih × jh, jh + h 2 .
Lemma 5.2. Up to the extraction of a subsequence, the functions Du h,∆t converge a.e. to Dũ in Q, and in L α (Q) for any α ∈ [1, 43 ).

Proof. Since ũ is the unique weak solution of (5.3), for every sequence of smooth functions μδ and smooth initial data μ0 δ converging to μ and to μ0 respectively, in the weak- * sense of measures, (which, for instance, can be constructed by convolution), the smooth solutions Ũδ satisfying

∂ t Ũδ -ν∆ Ũδ = μδ in Q, Ũδ (0, •) = μ0 δ in T 2
will converge to ũ, e.g. in L α (0, T ; W 1,α (T 2 )) for any α ∈ [1, 4 3 ), see e.g. [START_REF] Boccardo | Nonlinear elliptic and parabolic equations involving measure data[END_REF]. We now consider the finite difference approximation

U n+1 δ,i,j -U n δ,i,j ∆t -ν(∆ h U n+1 δ ) i,j = µ n δ,i,j , (5.4)
with (U 0 δ ) i,j = μ0 δ (x i,j ) and µ n δ,i,j = μδ (t n , x i,j ). Let σ be a positive real number: let T σ be the piecewise linear function defined on R by T σ (z) = max(-σ, min(z, σ)). Let the grid function e on Q h,∆t be given by e n i,j = u n i,j -U n δ,i,j . Define also, for 0 ≤ i, j < N h ,

I n j = i s.t. max[|e n i+1,j | , |e n i,j |] ≤ σ ; J n i = j s.t. max[|e n i,j+1 | , |e n i,j |] ≤ σ .
For any s : 0 < s < 1 we have

h 2 ∆t N T -1 n=0 i,j |D h e n+1 i,j | s ≤ h 2 ∆t N T -1 n=0 i,j |D + 1 e n+1 i,j | s + |D + 2 e n+1 i,j | s .
The first term can be estimated as follows

h 2 ∆t N T -1 n=0 i,j |D + 1 e n+1 i,j | s ≤ h 2 ∆t N T -1 n=0 j i∈I n+1 j (D + 1 e n+1 i,j D + 1 T σ (e n+1 ) i,j ) s 2 + h 2 ∆t N T -1 n=0 j i ∈I n+1 j |D + 1 e n+1 i,j | s ≤T 1-s 2   h 2 ∆t N T -1 n=0 i,j D + 1 e n+1 i,j D + 1 T σ (e n+1 ) i,j   s 2 +   h 2 ∆t N T -1 n=0 i,j |D + 1 e n+1 i,j |   s   h 2 ∆t N T -1 n=0 j i ∈I n+1 j 1    1-s ≤T 1-s 2   h 2 ∆t N T -1 n=0 i,j D + 1 e n+1 i,j D + 1 T σ (e n+1 ) i,j   s 2 +   h 2 ∆t N T -1 n=0 i,j |D + 1 e n+1 i,j |   s   h 2 ∆t N T -1 n=0 i,j |e n+1 i+1,j | + |e n+1 i,j | σ   1-s .
Similarly we estimate the term with D + 2 using the set J n i , and overall we deduce that

h 2 ∆t N T -1 n=0 i,j |D h e n+1 i,j | s ≤ c T 1-s 2   h 2 ∆t N T -1 n=0 i,j D h e n+1 i,j • D h T σ (e n+1 ) i,j   s 2 + c D h e s L 1 (Q h,∆t ) e 1-s L 1 (Q h,∆t ) σ -(1-s) (5.5)
for some constant c only depending on s. We estimate the first term from the discrete equation

ν h 2 ∆t N T -1 n=0 i,j D h e n+1 i,j • D h T σ (e n+1 ) i,j = h 2 ∆t N T -1 n=0 i,j (f n i,j -µ n δ,i,j )T σ (e n+1 i,j ) -h 2 ∆t N T -1 n=0 i,j e n+1 i,j -e n i,j ∆t T σ (e n+1 i,j )
which implies, using that (x -y)T σ (x) ≥ Θ σ (x) -Θ σ (y) for the nonnegative and convex function Θ(s) =

s 0 T σ (r)dr, ν h 2 ∆t N T -1 n=0 i,j D h e n+1 i,j • D h T σ (e n+1 ) i,j ≤ σ f -µ δ L 1 (Q h,∆t ) + e 0 L 1 (T 2 h ) .
Therefore, we deduce from (5.5)

h 2 ∆t N T -1 n=0 i,j |D h e n+1 i,j | s ≤ c T 1-s 2 σ s 2 f -µ δ L 1 (Q h,∆t ) + e 0 L 1 (T 2 h ) s 2 + c D h e s L 1 (Q h,∆t ) e 1-s L 1 (Q h,∆t ) σ -(1-s) .
Taking the minimum of the right hand side w.r.t. σ, and using the L 1 bounds for µ δ , µ 0δ and the data in (5.2), we see that

|D h e| s L 1 (Q h,∆t ) ≤ c D h e θs L 1 (Q h,∆t ) e θ(1-s) L 1 (Q h,∆t )
for some c and θ depending on s but not on h or δ. Recalling the definition of e, and the estimate on the discrete gradient, we have proved that

|D h u -D h U δ | s L 1 (Q h,∆t ) ≤ c u -U δ θ(1-s) L 1 (Q h,∆t ) .
Hence,

|D h u -D h U δ | s L 1 (Q h,∆t ) ≤ c u h,∆t -U δ,h,∆t θ(1-s) L 1 (Q h,∆t ) , (5.6)
where u h,∆t has been defined in Lemma 5.1 and U δ,h,∆t is the piecewise constant function that takes the value U n+1 δ,i,j in (t n , t n+1 )×(ih-h/2, ih+h/2)×(jh-h/2, jh+ h/2). Let us also define DU δ,h,∆t from the grid function U δ in a similar way as Du h,∆t : it takes the values

D + 1 U n+1 δ floor( i 2 ),ceil( j 2 ) , D + 2 U n+1 δ ceil( i 2 ),floor( j 2 ) in (t n , t n+1 ) × i h 2 , (i + 1) h 2 × j h 2 , (j + 1) h 2 .
Therefore, we see that

| Du h,∆t -Dũ| s L 1 (Q) ≤ | Du h,∆t -DU δ,h,∆t | s L 1 (Q) + | DU δ,h,∆t -D Ũδ | s L 1 (Q) + |D Ũδ -Dũ| s L 1 (Q) ≤c u h,∆t -U δ,h,∆t θ(1-s) L 1 (Q h,∆t ) + | DU δ,h,∆t -D Ũδ | s L 1 (Q) + |D Ũδ -Dũ| s L 1 (Q) ≤c u h,∆t -ũ θ(1-s) L 1 (Q h,∆t ) + ũ -Ũδ θ(1-s) L 1 (Q h,∆t ) + Ũδ -U δ,h,∆t θ(1-s) L 1 (Q h,∆t ) + | DU δ,h,∆t -D Ũδ | s L 1 (Q) + |D Ũδ -Dũ| s L 1 (Q)
where we have used (5.6) to obtain the third line. At fixed δ, since Ũδ is a smooth solution of the heat equation, the discrete approximation U δ,h,∆t converges to Ũδ in L 2 (Q) and DU δ,h,∆t converges to D Ũδ in L 2 (Q; R 2 ). Using also Lemma 5.1, we get that lim sup

h,∆t→0 |Du h,∆t -Dũ| s L 1 (Q) ≤ c ũ -Ũδ θ(1-s) L 1 (Q) + |D Ũδ -Dũ| s L 1 (Q) .
We conclude using the strong convergence of Ũδ to ũ in L α (0, T ; W 1,α (T 2 )) for any α ∈ [1, 4 3 ) (see e.g. [START_REF] Boccardo | Nonlinear elliptic and parabolic equations involving measure data[END_REF], as δ → 0. So

| Du h,∆t -Dũ| s L 1 (Q) ∆t,h→0
→ 0 which in particular implies that Du h,∆t converges to Dũ a.e. in Q and then, by Vitali's theorem, in L α (Q) for any α ∈ [1, 4 3 ). Remark 4. As a consequence of Lemma 5.2, for any ξ ∈ R 2 , there exists a subsequence of h and ∆t (not relabeled) such that the maps (t, x) → Du h,∆t (t, x+hξ) also converge to Dũ a.e. and in L α (Q) for any α ∈ [1, 4 3 ). Remark 5. Alternative strategies can be used to construct a function defined on Q from the grid function u. For example, we can define w h,∆t as the continuous and piecewise trilinear function on Q which takes the values u max(1,n) i,j at (t n , x i,j ) and which is trilinear in the rectangles of the time-space grid Q h,∆t . The advantage of taking w h,∆t instead of u h,∆t is that the former has weakly integrable partial derivatives with respect to the spatial variable. Therefore, we can use directly Dw h,∆t instead of having to define an independent approximation of Dũ such as Du h,∆t . It is then possible to prove the following lemma, which may replace both Lemmas 5.1 and 5.2:

Lemma 5.3. There exists a subsequence of h and ∆t (not relabeled) and a function ũ such that w h,∆t → ũ in L α (0, T ; W 1,α (T 2 )) for any α ∈ [1, 4 3 ) and in L β (Q) for all β ∈ [1, 2). In particular, w h,∆t → ũ and Dw h,∆t → Dũ in L 1 (Q) and almost everywhere in Q.

Proof. The strategy of proof is similar except that we may directly use the continuous version of the compactness lemma of Aubin-Simon, see [START_REF] Simon | Compact sets in the space L p (0, T ; B)[END_REF], for the function w h,∆t .

6. From the discrete to the continuous system. 6.1. A priori estimates and compactness. Let u h,∆t and Du h,∆t be the piecewise constant functions defined in Lemmas 5.1 and 5.2 respectively: up to the extraction of a subsequence, we can assume that u h,∆t → ũ in L β (Q) for all β ∈ [1, 2) and that Du h,∆t a.e. to Dũ in Q and in L α (Q) for any α ∈ [1, 4 3 ). Let m h,∆t be the piecewise constant function which takes the value m n i,j in (t n , t n+1 )× (ih -h/2, ih + h/2) × (jh -h/2, jh + h/2), and Dm h,∆t be the piecewise constant function from Q to R 2 which takes the value

D + 1 m n floor( i 2 ),ceil( j 2 ) , D + 2 m n ceil( i 2 ),floor( j 2 ) in (t n , t n+1 ) × i h 2 , (i + 1) h 2 × j h 2 , (j + 1) h 2 .
From Theorem 4.7, we may also assume that m h,∆t → m in L 1 (Q) and almost everywhere in Q, and that Dm h,∆t → D m weakly in L α (Q) for any α ∈ [1, 4 3 ). Moreover, for all η > 0, there exists a constant c η such that for all z ≥ 0, F (z) ≤ zF (z) η + c η . This fact and estimate (4.8) yield the equi-integrability of F (m h,∆t ). By Vitali's theorem, F (m h,∆t ) → F ( m) in L 1 (Q). From the observations above, the piecewise constant function which takes the value

u n+1 i,j -u n i,j ∆t -ν(∆ h u n+1 ) i,j -F (m n i,j )
in (t n , t n+1 )×(ih-h/2, ih+h/2)×(jh-h/2, jh+h/2) converges to ∂ ũ ∂t -ν∆ũ-F ( m) in the sense of distributions.

6.2. Stability of the discrete Bellman equation. We now pass to the limit in the discrete Bellman equation.

The main difficulty is to handle the nonlinear term g(x i,j , ∇ h u n+1 i,j ); here we wish to use the a.e. convergence of the gradients obtained in § 5. We adapt the method used for continuous problems in [START_REF] Boccardo | Existence de solutions non bornées pour certaines équations quasi-linéaires[END_REF]. Note that ∇ h u n+1 i,j is the value taken by the piecewise constant function with values in R

g h,∆t → g(x, D 1 ũ, D 1 ũ, D 2 ũ, D 2 ũ) = H(x, Dũ) a.e. in Q,
where g h,∆t is the piecewise constant function which take the value g(x i,j , ∇ h u n+1 i,j ) for (t, x) such that |x 1 -ih| < h/2, |x 2 -jh| < h/2, t n ≤ t < t n+1 .

Let now ϕ be a smooth function on T 2 such that ϕ ≥ 0, with ϕ(T ) = 0. We multiply the discrete Bellman equation by ϕ(t n+1 , x i,j ) and sum for all i, j and n = 0, . . . , N T -1. Since, by convexity, g(x, q) ≥ g(x, 0) + g q (x, 0) • q the regularity of g w.r.t. x and the L 1 -compactness of Du h,∆t allow us to apply Fatou's lemma obtaining

lim inf h→0 h 2 ∆t N T -1 n=0 i,j g(x i,j , ∇ h u n+1 i,j )ϕ(t n+1 , x i,j ) ≥ Q H(x, Dũ)ϕ dxdt .
Passing to the limit in the other terms of the equation, we deduce that

- Q ũ ϕ t dxdt+ν Q DũDϕ dxdt+ Q H(x, Dũ)ϕ dxdt ≤ T 2 u 0 ϕ(0) dx+ Q F ( m)ϕ dxdt .
We now wish to obtain the reverse inequality, which is the difficult part. We start by noticing that, since the monotonicity assumption implies g(x, q 1 , q 2 , q 3 , q 4 ) ≤ g(x, -q - 1 , q + 2 , -q - 3 , q + 4 ) from (3.1) and (3.2) and the fact that g(x, 0) is bounded, we know there exists λ > 0 such that g(x, q 1 , q 2 , q 3 , q 4 ) ≤ ν λ 1 + (q - 1 ) 2 + (q + 2 ) 2 + (q - 3 ) 2 + (q + 4 ) 2 . (

We multiply the discrete Bellman equation by e -λu n+1 i,j ϕ(t n+1 , x i,j ) and sum for all i, j and n = 0, . . . , N T -1. We obtain

h 2 ∆t N T -1 n=0 i,j u n+1 i,j -u n i,j ∆t e -λu n+1 i,j ϕ(t n+1 , x i,j ) + ν h 2 ∆t N T -1 n=0 i,j D h u n+1 i,j • D h (e -λu n+1 i,j ϕ(t n+1 , x i,j )) i,j + h 2 ∆t N T -1 n=0 i,j g(x i,j , ∇ h u n+1 i,j )e -λu n+1 i,j ϕ(t n+1 , x i,j ) = h 2 ∆t N T -1 n=0 i,j F (m n i,j )e -λu n+1 i,j ϕ(t n+1 , x i,j ). (6.2)
Since u is uniformly bounded below, the last term converges by dominated convergence, so

lim h→0 h 2 ∆t N T -1 n=0 i,j F (m n i,j )e -λu n+1 i,j ϕ(t n+1 , x i,j ) = Q F ( m)e -λũ ϕ dxdt. (6.3)
By convexity of s → e -λs and since φ(T, •) = 0, we have

h 2 ∆t N T -1 n=0 i,j u n+1 i,j -u n i,j ∆t e -λu n+1 i,j ϕ(t n+1 , x i,j ) ≤ 1 λ h 2 ∆t N T -1 n=0 i,j e -λu n i,j -e -λu n+1 i,j ∆t ϕ(t n+1 , x i,j ) = 1 λ h 2 ∆t N T -1 n=0 i,j e -λu n i,j ϕ(t n+1 , x i,j ) -ϕ(t n , x i,j ) ∆t + 1 λ h 2 i,j
e -λu 0 i,j ϕ(0, x i,j ) , and so, again by dominated convergence, lim sup

h→0 h 2 ∆t N T -1 n=0 i,j u n+1 i,j -u n i,j ∆t e -λu n+1 i,j ϕ(t n+1 , x i,j ) ≤ 1 λ Q e -λu ϕ t dxdt + 1 λ T 2 e -λu0 ϕ(0) dx . (6.4)
Let us deal now jointly with the second and third term in (6.2). First we split the energy term according to the sign of D + 1 (u) i,j (and D + 2 (u) i,j , respectively); indeed, we can write

D + 1 (u n+1 ) i,j D + 1 (e -λu n+1 ϕ n+1 ) i,j = (D + 1 (u n+1 ) i,j ) + e -λu n+1 i+1,j -e -λu n+1 i,j ϕ n+1 i+1,j -(D + 1 (u n+1 ) i,j ) -e -λu n+1 i+1,j -e -λu n+1 i,j ϕ n+1 i,j + (D + 1 (u n+1 ) i,j ) + (ϕ n+1 i+1,j -ϕ n+1 i,j )e -λu n+1 i,j -(D + 1 (u n+1 ) i,j ) -(ϕ n+1 i+1,j -ϕ n+1 i,j )e -λu n+1 i+1,j
and the same for the term with D + 2 . Reordering the indexes in the sum, this means that the D 1 part in the second order term can be read as

ν h 2 ∆t N T -1 n=0 i,j (D + 1 (u n+1 ) i-1,j ) + 1 h e -λu n+1 i,j -e -λu n+1 i-1,j ϕ n+1 i,j -ν h 2 ∆t N T -1 n=0 i,j (D + 1 (u n+1 ) i,j ) -1 h e -λu n+1 i+1,j -e -λu n+1 i,j ϕ n+1 i,j + ν h 2 ∆t N T -1 n=0 i,j (D + 1 (u n+1 ) i,j ) + D + 1 (ϕ n+1 ) i,j e -λu n+1 i,j -ν h 2 ∆t N T -1 n=0 i,j (D + 1 (u n+1 ) i,j ) -D + 1 (ϕ n+1 ) i,j e -λu n+1 i+1,j , which is equal to ν h 2 ∆t N T -1 n=0 i,j |(D + 1 (u n+1 ) i-1,j ) + | 2 e -λu n+1 i,j -e -λu n+1 i-1,j u n+1 i,j -u n+1 i-1,j ϕ n+1 i,j + ν h 2 ∆t N T -1 n=0 i,j |D + 1 (u n+1 ) i,j ) -| 2 e -λu n+1 i+1,j -e -λu n+1 i,j u n+1 i+1,j -u n+1 i,j ϕ n+1 i,j + ν h 2 ∆t N T -1 n=0 i,j (D + 1 (u n+1 ) i,j ) + e -λu n+1 i,j -(D + 1 (u n+1 ) i,j ) -e -λu n+1 i+1,j D + 1 (ϕ n+1 ) i,j .
We proceed similarly for the part with D 2 . Therefore,

ν h 2 ∆t N T -1 n=0 i,j D h u n+1 i,j • D h (e -λu n+1 i,j ϕ(t n+1 , x i,j )) i,j + h 2 ∆t N T -1 n=0 i,j g(x i,j , ∇ h u n+1 i,j )e -λu n+1 i,j ϕ(t n+1 , x i,j ) = ν h 2 ∆t N T -1 n=0 i,j (D + 1 (u n+1 ) i,j ) + e -λu n+1 i,j -(D + 1 (u n+1 ) i,j ) -e -λu n+1 i+1,j D + 1 (ϕ n+1 ) i,j + ν h 2 ∆t N T -1 n=0 i,j (D + 2 (u n+1 ) i,j ) + e -λu n+1 i,j -(D + 2 (u n+1 ) i,j ) -e -λu n+1 i,j+1 D + 2 (ϕ n+1 ) i,j + ν h 2 ∆t N T -1 n=0 i,j |(D + 1 (u n+1 ) i-1,j ) + | 2 e -λu n+1 i,j -e -λu n+1 i-1,j u n+1 i,j -u n+1 i-1,j ϕ n+1 i,j + ν h 2 ∆t N T -1 n=0 i,j |D + 1 (u n+1 ) i,j ) -| 2 e -λu n+1 i+1,j -e -λu n+1 i,j u n+1 i+1,j -u n+1 i,j ϕ n+1 i,j + ν h 2 ∆t N T -1 n=0 i,j |(D + 2 (u n+1 ) i,j-1 ) + | 2 e -λu n+1 i,j -e -λu n+1 i,j-1 u n+1 i,j -u n+1 i,j-1 ϕ n+1 i,j + ν h 2 ∆t N T -1 n=0 i,j |(D + 2 (u n+1 ) i,j ) -| 2 e -λu n+1 i,j+1 -e -λu n+1 i,j u n+1 i,j+1 -u n+1 i,j ϕ n+1 i,j + h 2 ∆t N T -1 n=0 i,j
g(x i,j , ∇ h u n+1 i,j )e -λu n+1 i,j ϕ n+1 i,j .

The first two terms in the right-hand side converge to ν T 0 T2 e -λũ Dũ • dφ by Lebesgue theorem, since Du h,∆t converges strongly in L 1 , ϕ is smooth and e -λu h is uniformly bounded and a.e. convergent. As far as the remaining terms are concerned, using that e -λs -e -λs s -s ≤ -λe -λ max(s,s ) , and due to (6.1), we observe that the last five terms under summation are bounded above, so that we can again apply Fatou's lemma, on account of the a.e. convergence of u h,∆t and Du h,∆t . Therefore, we conclude that lim sup for every smooth ϕ ≥ 0. In order to conclude, we need now to get rid of the exponential in the above inequality (6.6). To this purpose, we first observe that e -λũ ∈ L 2 (0, T ; H 1 (T 2 )) ∩ L ∞ (Q). (6.7) This can be easily proved obtaining an a priori estimate on e -λu h,∆t . Indeed, whenever u is a grid function which solves (5.1) for some data satisfying (5.2), we have

h 2 ∆t N T -1 n=0 i,j u n+1 i,j -u n i,j ∆t ψ(u n+1 i,j ) + ν h 2 ∆t N T -1 n=0 i,j D h u n+1 i,j • D h ψ(u n+1 i,j ) ≤ C ψ ∞
for any bounded real function ψ(r). In particular, if ψ is nondecreasing, this implies where C only depends on the L 1 -norm of the data. This is the desired a priori estimate; from which, using Fatou's lemma, we deduce

Q |Dũ| 2 ψ (ũ) dxdt ≤ C ψ ∞ .
On account of the fact that u is bounded below, we can take for example ψ(r) = 1 -e -µr to deduce that e -µũ ∈ L 2 (0, T ; H 1 (T 2 )) for any µ > 0. Thanks to (6.7), inequality (6.6) holds true not only for smooth functions ϕ but also for ϕ ∈ H 1 (Q) ∩ L ∞ , through a standard density argument. Moreover, there is no loss of generality in assuming that u 0 ∈ H 1 (T 2 ), so we extend ũ for negative t as identically equal to u 0 . Then, we choose ϕ(x, t) = ξ(t) 1 h e -λ(r-T k (r)) dr .

Moreover, 1 h t t-h e λT k (ũ)(x,s) ds converges to e λT k (u) in L 2 (0, T ; H 1 (T 2 )) and weak- * in L ∞ (Q), so we can pass to the limit as h → 0 in the remaining terms of (6.6). Finally, we obtain We conclude by letting k → ∞, thanks to the dominated convergence theorem:

- Q ũ ξ t dxdt- T 2 u 0 ξ(0) dx+ν Q DũDξ dxdt+ Q H(x, Dũ)ξ dxdt ≥ Q F ( m) ξ dxdt ,
for every ξ ≥ 0. Since the reverse inequality was already obtained previously, in the end, we proved that u solves

- Q ũ ξ t dxdt- T 2 u 0 ξ(0) dx+ν Q DũDξ dxdt+ Q H(x, Dũ)ξ dxdt = Q F ( m) ξ dxdt
for every ξ ∈ C 1 c ([0, T )), ξ ≥ 0, and therefore for every ξ. This concludes the proof that ũ is a weak solution to the limit equation.

6.3. Stability of the discrete Fokker-Planck equation. We now pass to the limit in the discrete Fokker-Planck equation.

By (4.6), the L 1 -compactness of m h,∆t and of Du h,∆t , we deduce the strong convergence in L 1 (Q) for the piecewise constant function which takes the value m n i,j ∇ q g(x i,j , [∇ h u n+1 ] i,j ) for (t, x) such that |x 1 -ih| < h/2, |x 2 -jh| < h/2, t n ≤ t < t n+1 . Moreover, by the consistency assumption we have,

h 2 ∆t N T -1 n=0 i,j m n i,j ∇ q g(x i,j , [∇ h u n+1 ] i,j ) • [∇ h ϕ] i,j → → Q ∇ q g(x, D 1 ũ, D 1 ũ, D 2 ũ, D 2 ũ) • (D 1 ϕ, D 1 ϕ, D 2 ϕ, D 2 ϕ)dxdt = Q m ∂H ∂p (x, Dũ) • Dϕ dxdt
Therefore, we can pass to the limit in the weak formulation and deduce that m is a weak solution of the Fokker-Planck equation.

We notice that the regularity m[H p (•, Dũ)Dũ -H(•, Dũ)] ∈ L 1 (Q) follows from inequality (4.10), by using Fatou's lemma. Moreover, we also find that m|H p (•, Dũ)| 2 ∈ L 1 (Q). The regularity ũ, m ∈ C 0 ([0, T ]; L 1 (T 2 )) follows from properties of weak solutions, see [START_REF] Porretta | Weak Solutions to Fokker-Planck Equations and Mean Field Games[END_REF]. Finally, this concludes the proof of Theorem 3.1.

  mT ) From these observations, (4.6), (4.7) and (4.8) follow from (4.10) and (3.1).

h→0 ν h 2 ∆tN T - 1 Q( 6 . 5 )F

 2165 n=0 i,j D h u n+1 i,j • D h (e -λu n+1 i,j ϕ(t n+1 , x i,j )) i,j + h 2 ∆t N T -1 n=0 i,j g(x i,j , ∇ h u n+1 i,j )e -λu n+1 i,j ϕ(t n+1 , x i,j ) ≤ ν DũDϕ e -λũ dxdt -ν λ Q |Dũ| 2 e -λũ ϕ dxdt + H(x, Dũ)e -λũ ϕ dxdt .Putting together (6.3)-(6.4)-(6.5), we deduce from (6.2) that ũ satisfies1 λ Q e -λũ ϕ t dxdt + 1 λ T 2 e -λu0 ϕ(0) dx + ν Q DũDϕ e -λũ dxdt -ν λ Q |Dũ| 2 e -λũ ϕ dxdt + Q H(x,Dũ)e -λũ ϕ dxdt ≥ Q ( m)e -λũ ϕ dxdt (6.6)

h 2 N T - 1

 21 n=0 i,j Ψ(u n+1 i,j ) -Ψ(u n i,j ) + ν h 2 ∆t N T -1 n=0 i,j D h u n+1 i,j • D h ψ(u n+1 i,j ) ≤ C ψ ∞ where Ψ(s) = s 0 ψ(r)dr. Thus, since |Ψ(s)| ≤ ψ ∞ |s|, one gets ν h 2 ∆t N T -1 n=0 i,j D h u n+1 i,j • D h ψ(u n+1 i,j ) ≤ C ψ ∞

ee

  λT k (ũ)(x,s) ds where T k (r) = min(r, k) and ξ ∈ C 1 c [0, T ). Using the monotone character of s → e -λs we have (see Lemma 2.3 in[START_REF] Blanchard | Stefan problems with nonlinear diffusion and convection[END_REF])lim sup h→0 1 λ Q e -λu ϕ t dxdt + 1 λ T 2 e -λu0 ϕ(0) dx≤ --λ(r-T k (r)) dr -

eeQF

  -λ(r-T k (r)) dr --λ(r-T k (r)) dr+ ν DũDξ e -λ(ũ-T k (ũ)) dxdt -ν λ {ũ>k} |Dũ| 2 e -λ(ũ-T k (ũ)) ξ dxdt + Q H(x, Dũ)e -λ(ũ-T k (ũ)) ξ dxdt ≥ Q ( m)e -λ(ũ-T k (ũ)) ξ dxdt.

  4 e 1 • Du h,∆t (• + h 2 e 1 ), e 1 • Du h,∆t (• -h 2 e 1 ), e 2 • Du h,∆t (• + h 2 e 2 ), e 2 • Du h,∆t (• -h 2 e 2 ) at (t, x) such that |x 1 -ih| < h/2, |x 2 -jh| < h/2, t n ≤ t < t n+1 . From the continuity of g, the consistency assumption and Remark 4,

il carattere C 1 rispetto a x dove viene realmente usato ? Idem per la g

Acknowledgements. The first author was partially funded by the ANR projects ANR-12-MONU-0013 and ANR-12-BS01-0008-01.