J Nghiem

E Giard

M Delmas

J B Rodriguez

P Christol

M Caes

H Martijn

E Costard

I Ribet- Mohamed

RADIOMETRIC CHARACTERIZATION OF TYPE-II INAS/GASB SUPERLATTICE (T2SL) MIDWAVE INFRARED PHOTODETECTORS AND FOCAL PLANE ARRAYS

I. INTRODUCTION

In recent years, Type-II InAs/GaSb superlattice (T2SL) has emerged as a new material technology suitable for high performance infrared (IR) detectors operating from Near InfraRed (NIR, 2-3µm) to Very Long Wavelength InfraRed (LWIR, λ > 15µm) wavelength domains. To compare their performances with well-established IR technologies such as MCT, InSb or QWIP cooled detectors, specific electrical and radiometric characterizations are needed: dark current, spectral response, quantum efficiency, temporal and spatial noises, stability…

In this paper, we first present quantum efficiency measurements performed on T2SL MWIR (3-5µm) photodiodes and on one focal plane array (320x256 pixels with 30µm pitch, realized in the scope of a french collaboration ). Different T2SL structures (InAs-rich versus GaSb-rich) with the same cutoff wavelength (λ c = 5µm at 80K) were studied. Results are analysed in term of carrier diffusion length in order to define the optimum thickness and type of doping of the absorbing zone. We then focus on the stability over time of a commercial T2SL FPA (320x256 pixels with 30µm pitch), measuring the commonly used residual fixed pattern noise (RFPN) figure of merit. Results are excellent, with a very stable behaviour over more than 3 weeks, and less than 10 flickering pixels, possibly giving access to long-term stability of IR absolute calibration.

II. InAs/GaSb T2SL QUANTUM STRUCTURE

A superlattice (SL) is a periodic stack of thin heterostructures, ie interfaces between two layers of different semiconductors. The InAs/GaSb heterointerface features a specific type-III band alignment, where the valence band of the GaSb layer is higher than the conduction band of the InAs layer, as shown in Fig 1a.

Each GaSb/InAs/GaSb block acts as a quantum well, thus this structure can be viewed as a system of coupled multi quantum wells. Electrons are then confined in minibands of energy instead of discrete levels. The bandgap of this periodic structure is defined by the energy difference between the first electron miniband C1 and the first heavy hole miniband V1 as illustrated in

III.QUANTUM EFFICIENCY MEASUREMENTS

A. Comparisons of several T2SL structure designs

Comparisons are first made with 3 monoelement structures that have the same active zone thickness of 1µm, as seen in Fig. 3. They were made by Molecular Beam Epitaxy (MBE) using GaSb as substrate and were designed in order to have the same wavelength cutoff (λ c =5µm). Sample A is GaSb-rich, with a period made of 10ML of InAs and 19 ML of GaSb (10/19). Sample B is symmetric (10/10). Sample C is InAs-rich (7/4). Quantum Efficiency (QE) measurements have been performed at 0V bias voltage on these 3 structures. First the relative photoresponse was measured using a FTIR spectrometer and then the absolute QE was calibrated using a SR200 blackbody [START_REF] Giard | Radiometric and noise characteristics of InAs-rich T2SL MWIR pin photodiodes[END_REF]. The QE results were compared in Fig. 3. At =4.5µm, better QE values were obtained by the symmetric and GaSb-rich SL structures while the absorption is higher in the InAs-rich design because of its large electronholes wave-function overlap. To explain why the QE of the InAs-rich structure is unexpectedly much lower than the other two, measurements at λ=4.5µm were analysed as a function of the bias voltage. Results are also shown in Fig 3 .  Increasing the bias voltage improves QE for the InAs-rich structure. It increases linearly from 6% to 17% by changing the bias voltage from 0V to -1V and then saturates. On the contrary, the QE of the GaSb-rich structures is independent of the bias voltage

The photocurrent is created in two steps: photon absorption (which creates the electron/hole pair) and collection of those carriers. The fact that increasing the bias voltage also improves the QE means that there is a problem of transport in the InAs-rich structure.

) ( * ) 1 ( 2 1 x x zce e e R QE        (2) ) ) cosh( ) sinh( ( * 1 ) 1 ( 2 max 2 max 2 2 1 2 e h x h h x e e p L L x x e L x x L e L L R QE                (3) p x h h x h h h h n e L L x e L x L L L R QE ) ) cosh( ) sinh( ( * 1 ) 1 ( 1 1 1 1 2 2              (4)
Where α is the absorption coefficient, set equal to 2200 cm -1 [START_REF] Tian | Influence of composition in InAs/GaSb type-II superlattices on their optical properties[END_REF]; L e and L h are the electron and hole diffusion lengths respectively, with L e set equal to 6µm [START_REF] Rehm | InAs/GaSb superlattices for advanced infrared focal plane arrays[END_REF]. R is the reflectance at the air/T2SL interface and set equal to 0.3. x 1 and x max -x 2 the thickness of the quasi neutral zone N and the quasi neutral zone P and depend on the bias voltage. The InAs-rich SL structure is n-type residual with hole minority carriers [START_REF] Cervera | Unambiguous determination of carrier concentration and mobility for InAs/GaSb superlattice photodiode optimization[END_REF]. We determined the hole diffusion length L h by fitting the equation with the measured QE, using L h as the only adjustable parameter. Fig 4 shows the QE simulation at V bias =0 and T =77K for three InAs-rich samples with three different active zone thicknesses equal to 500nm, 1µm and 4µm. The best agreement of the measured QE curves with Hovel equations was obtained with L h =80nm. This poor value in diffusion length does not allow the hole minority carrirs to reach the collection zone and then penalizes the QE.

B. Conclusions

In the InAs-rich structure, the active zone features a residual N-type doping, meaning that the minority carriers are holes. As holes have a short diffusion length, the collection process is not optimal. One solution would be to change the minority carriers into electrons, which have a better diffusion length [START_REF] Cervera | Unambiguous determination of carrier concentration and mobility for InAs/GaSb superlattice photodiode optimization[END_REF]. This can be done by slightly P-doping the active zone, leading to a P + P -N junction as used in [START_REF] Hoffman | Beryllium compensation doping of InAs/GaSb infrared superlattice photodiodes[END_REF] to optimize LWIR structures. Fig 5 displays the QE comparison between intentionally P doped and non intentionally doped (nid) N-type residual InAs-rich photodiodes (7/4, 4µm active zone thickness, front side illuminated). We note that the measured QE is enhanced by up to 10 times by P-doping the InAs-rich T2SL (40% QE efficiency at λ=4.5µm for the P doped structure). To complete the comparison, QE measurements as a function of bias voltage are also plotted. In the P doped structure, we can see that the QE does not vary with the bias voltage, staying at 40% for λ=4.5µm, highlighting an optimised collection of electron minority carriers in p-doped InAs-rich T2SL structure. Another solution to improve the QE is to ensure that the holes are created closer to the P/N junction, by backside illuminating the T2SL just like the FPA configuration. Fig 6 represents the QE measurement of an InAs-rich nid MWIR FPA (320x256) made in France [START_REF] Giard | Quantum efficiency investigations of type-II InAs/GaSb midwave infrared superlattice photodetectors[END_REF]. The 42% average QE is due to the fact that this structure is double-pass and that the depletion zone is now betterr overlapping with the absorption zone, which allows more hole minority carriers to be collected.

Fig 6 Left : QE measurements on InAs-rich MWIR FPA (320x256) at 77K on 600 pixels at the center; right : schematic view of the back-side illuminated InAs-rich photodiode (the p-n junction is now at the p/i interface, allowing a good overlap between the absorption and the collection zones).

IV. STABILITY OVER TIME AND CORRECTABILITY

Regardless of the applications, uniformity and stability over time of FPAs have been key priorities and recent works on flickering pixels [START_REF] Brunner | Improvement of RTS Noise in HgCdTe MWIR Detectors[END_REF] prove that they are still relevant. The usual approach to evaluate them is through estimating the Residual Fixed Pattern Noise (RFPN) and the number of flickering pixels.

A. Residual Fixed Pattern Noise

Spatial noise is due to various physical effects that can appear in the photodiode itself, in the readout circuit or during the technological steps of hybridization. It results in the fact that pixels' response is not absolutely identical in an FPA. By using a gain/offset correction table, it is possible to damp the impact of those non-uniformities.

In order to obtain the gain/offset table, we used a Two Point Correction (TPC). It consists in the use of a blackbody at two fixed fluxes, or rather at two temperatures T bb with radiances L bb given by Planck's law for a blackbodies:

   d e hc T L kTbb hc bb Tbb    1 1 2 ) ( 5 2 (5)
where h is the Planck constant,k is the Boltzmann constant, c is the speed of light and λ is the wavelength. Gain and offset for each pixel are calculated such that all the pixels deliver the same signal when exposed to blackbodies at T1 and T2. As the pixels' behavior is not rigorously linear, the table calculated will not be absolutely correct at temperature between T1 and T2. This explains the W shape seen on measurements of RFPN, which is the remaining spatial noise after correction. RFPN has been measured on commercial IRNova T2SL MWIR FPA (320x256 FPAs with 30μm pixel pitch). Measurements have been done using the following routine:

 Cooling down the FPA in the morning.  One measurement when the FPA reaches its operating temperature (80K).

 One measurement at the end of the day.  The cooling machine is turned off, bringing the FPA back to room temperature.

ICSO 2016

Biarritz, France International Conference on Space Optics 18 -21 October 2016

This makes 2 measurements a day, except for the first day when more measurements were done.

For each measurement, gain and offset corrections were applied on each pixel, using the G/O tables determined on day D0. We consider the correction is still valid as long as RFPN is lower than the temporal noise (TN) measured at a given day.

In the case of the IRNova FPA, results in Fig 7 show RFPN/TN for each measurement. We can see that even after 3 weeks, the correction can be deemed as valid, which means that the stability over time is excellent. The number of flickering pixels is also very low (typically 10 pixels across the full FPA). This result means that a calibration can be made and does not need to be refreshed for a long time, which ultimately increases the IRFPA uptime (as the IR device has to be disabled during calibration phases),. In some cases, the IR device could even be designed without on-board calibration device, which greatly impacts the Size Weight and Power (SWaP) issue.

Fig 7 RFPN/rms using the same correction on different days

B. Discussion

Results in terms of RFPN and number of flickering pixels seem very promising for T2SL technology, potentially ensuring long-term stability of infrared instruments absolute calibration.

In a more general manner, one may wonder whether RFPN is the most relevant figure of merit to faithfully depict stability over time. Indeed, RFPN highlights the nonlinearity of pixels. It may not be a major concern as long as the behavior is stable, thus correctable. This preliminary work shows that only using RFPN standards will not be enough to state on the "correctability" (or ability to be calibrated on the long-term) of the FPA. For example, flickering pixels (which may bias the "standard" pixel response used for corrections), may be overlooked as the RFPN cannot discriminate them (and apply proper correction to them). We may reach a scenario where RFPN is very low (if no flickering occurred during the calibration) and the resulting image however cannot be considered as stable. Complementary measurements are ongoing at ONERA, with the objective to define a measurement protocol to quantify the "correctability" of FPAs and to compare the merits of different competing technologies. The electro-optic characterizations realized so far by different laboratories [START_REF] Giard | Quantum efficiency investigations of type-II InAs/GaSb midwave infrared superlattice photodetectors[END_REF][10] [11] [12][13] [START_REF] Höglund | Minority carrier lifetime and photoluminescence studies of antimony-based superlattices[END_REF] allow us to conclude that MWIR photodiodes and FPA exhibit a dark current higher than initially expected, typically 1000 times higher than the MCT state-of-the-art given by the Rule07 [9], making High Operating Temperature (HOT) MWIR out of reach for the time being. This high dark current has been attributed to low minority carriers lifetime (in the 80-100 ns range for T2SL in the MWIR domain), which is also responsible for the 80nm holes diffusion length reported in this paper. However, we have shown that a proper design of the active zone mitigates the impact of this low lifetime on QE values. The quantum efficiency of photodiodes is in the range of 50%, which is a very satisfying value, even if not as high as the 80% QE published in MCT technology. T2SL also proved to have no excess noise and flat angular response as well [START_REF] Giard | Radiometric and noise characteristics of InAs-rich T2SL MWIR pin photodiodes[END_REF]. We also pointed out that the long-term stability of T2SL technology is excellent, such that T2SL could get back into the race for high performance infrared detectors in the MWIR domain.

However, we believe that it is interesting for T2SL to focus on the LWIR (8-12µm) and VLWIR (>12µm) spectral ranges. Indeed, there is a need for a technology that would offer good E-O performances (especially QE) with an excellent stability, few flickering pixels, and available in large format FPA. We believe that T2SL can address this need, since very promising results have been already reported with a LWIR megapixel T2SL FPA (λ c =11µm@77K with 80% QE) [START_REF] Manurkar | High performance long wavelength infrared mega-pixel focal plane array based on type-II superlattices[END_REF] and a LWIR/VLWIR bispectral 320x256 FPA (λ c =11µm@77K with QE in the range of 40% in both bands) [START_REF] Huang | Type-II superlattice dual-band LWIR imager with M-barrier and Fabry-Perot resonance[END_REF]. LWIR T2SL photodiodes are currently being processed at IES and should be available for E-O characterization at ONERA very soon.

V.CONCLUSION

In conclusion, we studied different T2SL MWIR photodiodes (from InAs-rich to GaSb-rich SL structures) in terms of quantum efficiency. Measurements pointed out a carrier collection issue as higher bias voltage improves the quantum efficiency in InAs-rich structure. Matching the results in the InAs-rich with Hovel equations, we found that the collection zone is smaller than the active zone because of low carrier diffusion lengths. Solutions have been proposed: 1) changing the minority carrier from hole to electrons which have better diffusion length, or 2) switching to backside illumination in order to match the active zone with the collection zone. Experiments on an IRNova MWIR commercial T2SL FPA show that the stability over time is excellent based on RFPN measurements. It potentially ensures long-term stability of infrared instruments absolute calibration. This demonstrates that T2SL technology is competitive for high performance MWIR photodetection. Finally, "Correctability" has been introduced as an alternative to RFPN measurements, aiming to depict the stability of FPAs over extended periods of time.

  Fig 1, and depends on the thicknesses of the InAs and GaSb layers.

Fig 1 Fig 2

 12 Fig 1 Schematic InAs/GaSb SL structure displaying the fundamental electron C1 and heavy hole V1 minibands. The absorption of infrared photons occurs between these two minibands.

Fig 3

 3 Fig 3 Right :Comparison of quantum efficiency of 3 different T2SL photodiodes at 77K. Left : QE as a function of the wavelength, for Vbias=0V ; right : QE as a function of the bias, for λ=4.5µm.

Fig 4

 4 Fig 4 Left : comparison between QE measurements and predictions from Hovel equations for three InAs-rich photodiodes with different active zone thicknesses (500 µm, 1 µm and 4 µm), at 77K and 0V bias voltage; Right : schematic view of the front-side illuminated InAs-rich photodiode.

Fig 5

 5 Fig 5 Left : QE measurements on InAs-rich T2SL photodiodes with non intentionally doped or P doped active zone at 77K and V bias =0V ; Right : schematic view of the front-side illuminated InAs-rich photodiode (the p-n junction is now at the i/n interface, allowing a good overlap between the absorption and the collection zones).

  OF THE T2SL EMERGING TECHNOLOGY TO ADDRESS THE MWIR SPECTRAL DOMAINT2SL technology has initially drawn attention in the MWIR domain since it could potentially combine the advantages of MCT (high operating temperature) and InSb (excellent uniformity/stability) technologies, and thus potentially fulfill the SWaP (size, weight and power) requirements.

VI. ACKNOWLEDGEMENTS

The authors acknowledge the financial support of the French procurement agency (DGA) and the LabEx FOCUS.


Vanel Steve

Siyou Fotso

Engelbert Mephu Nguifo

Philippe Vaslin

Engelbert Mephu Nguifo

Parameter Free Piecewise Dynamic Time Warping for time series classification

Several improvements have been done in time series classification over the last decade. One of the best solutions is to use the Nearest Neighbour algorithm with Dynamic Time Warping(DTW), as the distance measure. Computing DTW is relatively expensive especially with very large time series. Piecewise Dynamic Time Warping (PDTW) is an efficient variant which consists of segmenting time series into fixedlength segments. However, the choice of the optimal size (or number) of segments remains a difficult challenge for end users. The Brute-force solution, a naive solution, repeats the classification with each segment size, and selects the one with the best accuracy. This solution is not appropriated especially when dealing with massive and large time series data. In this work, we propose a parameter free approach for PDTW, that finds the size (or number) of segments to be used with the Nearest Neighbour algorithm. Our approach is a heuristic that is parameter free since it does not require any domain specific tuning. Several properties of our heuristic are studied, and an extensive experimental comparison demonstrates its efficiency and effectiveness, in terms of accuracy and runtime.

Introduction

Time series are ubiquitous in sciences as for example in economics [START_REF] Hamilton | A new approach to the economic analysis of nonstationary time series and the business cycle[END_REF], in medicine [START_REF] Huang | Ecg frame classification using dynamic time warping[END_REF], in finance [START_REF] Marszalek | Modeling and forecasting financial time series with ordered fuzzy candlesticks[END_REF] or in computer science [START_REF] Myers | Performance tradeoffs in dynamic time warping algorithms for isolated word recognition[END_REF]. An important task is time series comparison that can be done in two main ways. Either the comparison method considers that there is no time distortion as in Euclidian distance (ED), or it considers that some small time distortions exist between time axis of time series as in Dynamic Time Warping alignment algorithm (DTW) [START_REF] Zhang | Dynamic time warping under pointwise shape context[END_REF]. Since time distortion often exists between time series, DTW often has better results than the ED [START_REF] Chen | The ucr time series classification archive[END_REF]. An exhaustive comparison of time series algorithms [START_REF] Bagnall | Time Se-ries Classification Website[END_REF] shows that DTW is among the efficient techniques to be used. However, DTW has two major drawbacks: the comparison of two time series under DTW is time-consuming [START_REF] Rakthanmanon | Searching and mining trillions of time series subsequences under dynamic time warping[END_REF] and DTW sometimes produces pathological alignments [START_REF] Keogh | Derivative dynamic time warping[END_REF]. A pathological alignment occurs when, during the comparison of two time series X and Y , one datapoint of the time series X is compared to a large subsequence of datapoints of Y . A pathological alignment causes a wrong comparison.

Three categories of methods are used to avoid pathological alignments with DTW:

• The first one adds constraints to DTW [START_REF] Ratanamahatana | Making Time-series Classification More Accurate Using Learned Constraints[END_REF], [START_REF] Yu | Dynamic time warping constraint learning for large margin nearest neighbor classification[END_REF], [START_REF] Candan | sdtw: computing dtw distances using locally relevant constraints based on salient feature alignments[END_REF], [START_REF] Sakoe | Dynamic programming algorithm optimization for spoken word recognition[END_REF], [START_REF] Jeong | Weighted dynamic time warping for time series classification[END_REF]. The main idea here is to limit the length of the subsequence of a time series that can be compared to a single datapoint of another time series.

• The second one suggests to skip datapoints that produce pathological alignment during the comparison of two time series [START_REF] Longin | Elastic partial matching of time series[END_REF], [START_REF] Itakura | Minimum prediction residual principle applied to speech recognition[END_REF], [START_REF] Myers | Performance tradeoffs in dynamic time warping algorithms for isolated word recognition[END_REF].

• The third one proposes to replace the datapoints of time series by a high level abstraction that captures the local behavior of those time series. A highlevel abstraction can be a histogram of values that captures the repartition of time series datapoints in space [START_REF] Zhang | Dynamic time warping under pointwise shape context[END_REF], or a feature that captures the local properties of time series, such as the trend with Derivative DTW (DDTW) [START_REF] Keogh | Derivative dynamic time warping[END_REF] or the mean with Piecewise DTW (PDTW) [START_REF] Keogh | Scaling up dynamic time warping for datamining applications[END_REF].

PDTW has been introduced with the aim to speed up the computation of DTW, which depends on the length of the time series. PDTW suggests to use a compact abstraction of time series instead of the raw data. Indeed, PDTW proposes to split a time series into consecutive fixed-length segments and to compute the mean of each segment. Then, the mean is used instead of the data points in the segment to compare the time series.

In practice, a straightforward way to use PDTW is the brute-force approach that consists in exploring all the possible values for the number of segments. However, this is not feasible with long time series data. So, the question is how to automatically fix this parameter without a considerable decrease of classification accuracy ?

In this paper, we propose a parameter free heuristic to align piecewise aggregate time series with DTW that approximates the optimal value of the number of segments to be considered during the alignment. In this heuristic, the number of segments is chosen based on the quality of the alignment, which is evaluated by the classification error on the training set. The best classification algorithm to use for this purpose is one Nearest Neighbor (1NN) that is combined with PDTW. In this case 1NN is the best because its classification error directly depends on the alignment of time series, since it has no other parameters [START_REF] Wang | Experimental comparison of representation methods and distance measures for time series data[END_REF].

2 Background and related work

Definition 2.1. A time series X = x 1 , • • • , x n
is a sequence of numerical values representing the evolution of a specific quantity during the time. x n is the most recent value.

2.1 Dynamic Time Warping algorithm. DTW [START_REF] Keogh | Exact indexing of dynamic time warping[END_REF] is a time series alignment algorithm that performs a non-linear alignment while minimizing the distance between two time series. To align two time series :

X = x 1 , x 2 , • • • , x n ; Y = y 1 , y 2 , • • • , y m .
the algorithm constructs an n × m matrix where the cell (i, j) of the matrix corresponds to the squared distance (x i -y j ) 2 that is the alignment between x i and y j . To find the best alignment between two time series, it constructs the path that minimizes the sum of squared distances. This path, noted W = w 1 , w 2 , . . . , w k , . . . , w K , should respect the following constraints:

• Boundary constraint: w 1 = (1, 1) and w K = (n, m)

• Monotonicity constraint: Given w k = (i, j), w k+1 = (i , j ) then i ≤ i and j ≤ j

• Continuity constraint: Given w k = (i, j), w k+1 = (i , j ) then i ≤ i + 1 and j ≤ j + 1

The warping path is computed by using an algorithm based on the dynamic programming paradigm that solves the following recurrence:

γ(i, j) = d(x i , y j )+min{γ(i-1, j-1), γ(i-1, j), γ(i, j-1)},
where d(x i , y j ) is the squared distance contained in the cell (i, j) and γ(i, j) is the cumulative distance at the position (i, j) that is computed by the sum of the squared distance at the position (i, j) and the minimal cumulative distance of its three adjacent cells. Definition 2.2. A segment X i of length l of the time series X of length n (l < n) is a sequence constituted by l variables of X starting at the position i and ending at the position i + l -1. We have: X i = x i , x i+1 , ..., x i+l-1 Definition 2.3. The arithmetic average of the data points of a segment X i of length l is noted Xi and is defined by:

Xi = 1 l l-1 j=0
x i+j Definition 2.4. Let T be the set of time series. The Piecewise Aggregate Approximation (PAA) is defined as follows:

P AA : T × N * → T (X, N ) → P AA(X, N ) = X1 , • • • , XN if N < |X| X otherwise
Piecewise Dynamic Time Warping Algorithm (PDTW) [START_REF] Keogh | Scaling up dynamic time warping for datamining applications[END_REF] is the DTW algorithm applied on Piecewise Aggregate time series [START_REF] Keogh | Dimensionality Reduction for Fast Similarity Search in Large Time Series Databases[END_REF]. Let N ∈ N * , X and Y be two time series.

P DT W (X, Y, N ) = DT W (P AA(X, N ), P AA(Y, N )).
The number of segments N that one considers greatly influences the quality of the alignment of the time series. However, PDTW does not give any information on the way to choose it. To do so, [START_REF] Chu | Iterative deepening dynamic time warping for time series[END_REF] proposes the Iterative Deepening Dynamic Time Warping Algorithm (IDDTW).

Iterative

Deepening Dynamic Time Warping. IDDTW only considers values for the number of segments that are power of 2 and for each value, computes an error distribution by comparing PDTW with the standard DTW at each level of compression. It takes as input: the query Q, the dataset D, the user's confidence (or tolerance for false dismissals) user conf , and the set of standard deviations StdDev obtained from the error distribution.

• The algorithm starts with applying the classic DTW to the first K candidates from the dataset. The results of the best matches to the query are contained in R, with |R| = K. The best so f ar is determined from argmaxR.

• Both the query Q and each subsequent candidate C are approximated using PAA representations with N segments to determine the corresponding PDTW.

• A test is performed to determine whether the candidate C can be pruned off or not. If the result of the test is found to have a probability that it could be a better match than the current best so f ar, a higher resolution of the approximation is required. Then each segment of the candidate is split into two segments to obtain a new candidate.

• The process of approximating Q and C to determine the PDTW should be reapplied and the test is repeated for all levels of approximations until they fail the test or their true distance DTW is determined.

Doing so, IDDTW finds the number of segments that best approximates DTW and speeds up its computation. However, the goal of IDDTW is not the same as ours, which is to find the number of segments that best aligns the time series and speeds-up the computation of DTW. Actually, IDDTW has three main drawbacks:

• It only considers the numbers of segments for PDTW that are power of 2;

• It requires a user-specified tolerance for false dismissals that influences the quality of the approximation, but the algorithm does not give any indication on how to choose the tolerance;

• It considers DTW as a reference while looking for the number of segments that best aligns the time series. However, because of pathological alignments, DTW sometimes fails to align time series properly.

In this paper, we propose a heuristic named parameter Free piecewise DTW (FDTW) that deals with all the drawbacks of IDDTW: it considers all the possible values for the number of segments, it is parameter-free and it finds a number of segments for PDTW based on the quality of the time series alignment namely the classification error. The next section presents a definition of our heuristic.

3 Heuristic search of the number of segments 3.1 Problem definition. Let D = {d i } be a set of datasets composed of time series. We note |d i | the number of time series of the dataset d i .

Let X ∈ d i be a time series of the dataset d i ; we note |X| = n the length of the time series X. For simplicity of notation we suppose that all the time series of d i have the same length.

Definition 3.1. 1N N DT W : D → [0, 1] d i → 1N N DT W (d i ) 1N N DT W (d i ) is
(d i , N ) → 1N N P DT W (d i , N ) = = 1N N DT W • P AAset(d i , N ) 1N N P DT W (d i , N
) is the classification error of 1-NN with PDTW using N segments on d i .

Our goal is to find the number of segments that allows P DT W to best align time series. P DT W gives a good alignment when its classification error with 1N N is low [START_REF] Rakthanmanon | Searching and mining trillions of time series subsequences under dynamic time warping[END_REF]. Our problem is then to find the number of segments N that minimizes 1N N P DT W (d i , N ).

Formaly, given a dataset d i , we look for the number of segments N ∈ {1 . . . n} such that

1N N P DT W (d i , N ) = min 1≤α≤n {1N N P DT W (d i , α)}.
3.2 Brute-force search. The simplest way to find the value for the number of segments that minimized the classification error is to test all the possible values. Obviously, this method is time consuming as we have to test n values to find the one that has the minimal classification error. The time complexity of this process is :

O(( |trainingset| 2 ) 2 × N ∈C N 2 ), |C| = n,
where C is the set of values for the number of segments.

To reduce the time of the search, the heuristic proposes to look for the number of segments with the minimal classification error without testing all the possible values.

3.3

Parameter free heuristic. The idea of our heuristic is the following:

1. We choose N c candidates distributed in the space of possible values to ensure that we are going to have small, medium and large values as candidates. The candidates values are: n, n-n Nc , n-2× n Nc , ..., n-N c × n Nc . For instance, if the length of time series is n = 12 and the number of candidates is N c = 4, we are going to select the candidates 12, 9, 6, 3.

1, 2, [3], 4, 5, [6], 7, 8, [START_REF] Itakura | Minimum prediction residual principle applied to speech recognition[END_REF], 10, 11, [START_REF] Keogh | Exact indexing of dynamic time warping[END_REF] 2. We evaluate the classification error with 1N N P DT W for each candidate that we have previously chosen and we select the candidate that has the minimal classification error: it is the best candidate. In our example, we may suppose that we get the minimal value with the candidate 6 it is thus the best candidate at this step.

1, 2, 3, 4, 5, [6], 7, 8, 9, 10, 11, 12

3. We respectively look between the predecessor (i.e., 3 here) and successor (i.e., 9 here) of the best candidate for a number of segments with a lower classification error. This number of segments corresponds to a local minimum. In our example, we are going to test the values 4, 5, 7 and 8 to see if there is a local minimum.

4. We restart at step one, while choosing differents candidates during each iteration to ensure that we return a good local minimum. We fix the number of iterations to log(n) . At each iteration the first candidate is n -(number of iteration -1).

In short, in the worst case, we test the N c first candidates to find the best one. Then, we test 2n Nc other candidates to find the local minimum. We finally perform nb(N c ) = N c + 2n

Nc tests. The number of tests that we have to perform is a function of the number of candidates. How many candidates should we consider to reduce the number of tests? The first derivative of the function nb vanishes when N c = √ 2n and the second derivative is positive so the minimal number of tests is done when the number of candidates N c = √ 2n. Algorithm 1 presents the details of the heuristic.

Time complexity: We use the training set to find the number of segments that should be considered with PDTW. To do so, we applied 1N N on the training set that costs A heuristic does not always give the optimal value. To ensure that it gives a result not far from the optimal value, one approach is to guarantee that the result of the heuristic always lies in an interval with respect to the optimal value [START_REF] Ibarra | Fast approximation algorithms for the knapsack and sum of subset problems[END_REF].

O(( |trainingset| 2 ) 2 × N ∈C N 2 ), |C| = log(n)-1 i=0 8 √ n -i.
In our case, we are looking for the number of segments that allows a good alignment of time series. The alignment is good when the classification error with 1NN is minimal or when the accuracy is maximal.

Let d i be a dataset: acc max(di) = 1min 1≤α≤n {1N N P DT W (di, α)} is the maximal accuracy for the dataset d i , acc DT W = 1 -1N N DT W (d i ) is the accuracy with d i and 1NNDTW and acc F DT W = 1 -F DT W (d i ) is the accuracy of our heuristic.

To ensure the quality of our heuristic FDTW, Proposition 3.1 assume that 1N N DT W is better than the baseline classifier Zero Rule. Zero rule classifier is a simple classifier that predicts the majority class of test data (if nominal) or average value (if numeric). Zero rule is often used as baseline classifier [START_REF] Cuřín | Meeting state recognition from visual and aural labels[END_REF]. The minimal value of the accuracy of Zero rule is 1 c where c is the number of classes of the dataset. Proposition 3.1. For a given dataset

d i that has c i classes, c i ∈ N * , if acc DT W ≥ 1 ci then 1 ci × acc max ≤ acc F DT W ≤ acc max Proof. By definition, acc F DT W ≤ acc max We look for k ∈ N such that 1 k × acc max ≤ acc F DT W 1 k × acc max ≤ acc F DT W i.e. acc max acc F DT W ≤ k or acc max acc F DT W ≤ 1 acc F DT W because acc max ≤ 1 and 1 acc F DT W ≤ 1 acc DT W because acc DT W ≤ acc F DT W 1 acc DT W ≤ c i because 1 c i ≤ acc DT W

by hypotesis

We take k = c i 4 Experiments and discussion 4.1 Datasets. The performance of FDTW has been tested on 45 datasets of the UCR time series datamining archive [START_REF] Chen | The ucr time series classification archive[END_REF], which provides a large collection of datasets that cover various domains (Table 1). Each dataset is divided into a training set and a testing set. The 45 datasets possess between 2 and 50 classes, the length of the time series varies from 24 to 1882, the training sets contain between 20 and 1000 time series and the testing sets contain between 28 and 6164 time series. An implementation of BF, IDDTW and FDTW is available online [START_REF] Siyou Fotso | Source code of FDTW[END_REF] 4.2 Results. Firstly, to evaluate the quality of our heuristic FDTW, we compared its classification errors with that of IDDTW (Figure 4) and the minimal one (Figure 3). The classification error was calculated based on the holdout model evaluation and the minimal one was find by applying Brute-force search (BF) on both training set and testing set. FDTW and IDDTW used the training set to find the number of segments N using 3 fold cross validation. IDDTW tested all the values of N that were equal to a power of two and kept the one that had a minimum classification error. We also compared FDTW to BF and IDDTW in terms of number of tested values (Figure 1), running time (Figure 2) and compression ratio.

Then, we compared FDTW to other classification methods reported in the literature. The comparison was based on the classification error calculated using the hold out evaluation model. The smallest classification error reported on each dataset and the 1NN classification error of Euclidean distance, DTW without a warping window and DTW with best warping window have been published by previous researchers [START_REF] Chen | The ucr time series classification archive[END_REF] 1) clearly shows that the number of candidates in BF is considerably reduced in FDTW by a factor at least greater than 2.5. This number is exponentially correlated to the time series length, for example FDTW tested 0.08% less candidates than BF with the dataset ItalyPowerDemand that has the shortest time series length of our sample (24 data points) and 76% less candidates than BF with dataset Inlineskate that has the longest time series of our sample (1882 data points). Actually, the number of candidates to be tested ranges from 1 to n, n being the length of time series. This demonstrates an advantage of FDTW in terms of space exploration and thus indirectly in terms of mem-ory usage and execution time. However, FDTW tested more candidates than IDDTW, which tested in average 96% less candidates than Brute-force search (Figure 1). Generally, FDTW is 8 times faster than Brute-force search with an average execution time of 176 minutes against 1386 minutes for Brute-force search. IDDTW is 7 times faster than FDTW and remains the fastest with an average execution time of 24 minutes. The execution time increases with the length of the time series (Figure 2). The increase of Brute-force search execution time is faster than that of FDTW and IDDTW. This can be seen on the dataset Lightning-2 whose time series have a length equal to 637 data points.

As regard in the compression ratio, the heuristic uses a compact representation for time series whose length contains in average 44% data points less than the initial time series against 63% for IDDTW.

The experiments show that IDDTW is faster and test fewer candidates. However, FDTW have better performance. Actually, FDTW resulted in a lower classification error than IDDTW on 22 datasets and the same classification error than IDDTW on 8 datasets (Figure 4). They also show that the classification error of Brute-force search (BF) is smaller than the smallest classification error reported in the literature on six datasets and is equal to the smallest classification error reported on four datasets. Moreover, In average, BF is better than the other algorithms of Table 4.3 with a classification error of 0.175. In other words, it is a good strategy to piecewise aggregate time series before classifying them if we know a good number of segment to use.

Our heuristic FDTW managed to find the minimum error for 9 datasets (Coffee, ECGFiveDays, Gun-point, ItalyPowerDemand, OliveOil, Plane, Synthetic control, Trace, Two patterns). It also outperforms the smallest classification error reported in the literature on dataset CBF (N • 5).

The methods of the literature are

Conclusion

Our problem was to choose a good number of segments for Piecewise Dynamic Time Warping. To answer this question, we proposed a heuristic approach called Parameter Free Piecewise Dynamic Time Warping (FDTW) that proposes an approximation of the best number of segments to be used during times series classification based on DTW. FDTW has been experimented on 45 data sets on a classification task. In average, it returned a classification error lower than the one of ID-DTW. Our approach is a heuristic that is parameter free since it does not require any domain specific tuning. This work allows to reduce the storage space and the processing time of time series classification without decreasing the quality of the alignment. As a perspective, we plan to use piecewise aggregate time series with other variants of DTW to improve the classification. Using the same strategy presented in FDTW, we plan to find the number of segments to be considered for symbolic representations of time series like SAX [START_REF] Lin | A symbolic representation of time series, with implications for streaming algorithms[END_REF], ESAX [START_REF] Lkhagva | Extended sax: Extension of symbolic aggregate approximation for financial time series data representation[END_REF], SAX-TD [START_REF] Sun | An improvement of symbolic aggregate approximation distance measure for time series[END_REF].

Results reported in [3][1]

Our experiments

  the classification error of one Nearest Neighbour with Dynamic Time Warping on the dataset d i . Definition 3.2. Let d ⊆ T be a subset of time series, N ∈ N * , P AAset(d, N ) = {P AA(X, N ), ∀X ∈ d} Definition 3.3. 1N N P DT W : D × {1 . . . n} → [0, 1]

where ( |trainingset| 2 ) 2 Lemma 3 . 1 .

 2231 comes from 1N N algorithm and N ∈C N 2 comes from P DT W with N being one value of the number of segments, and C being the set of values for the number of segments. At each iteration, the heuristic tests nb( √ 2n) = 8 √ n number of segments.We have log(n) iterations so |C| = For a given a datasetd i F DT W (d i ) ≤ 1N N DT W (d i ).The quality of the alignment of our heuristic is better than that of DTW.Proof. 1N N DT W (d i ) = 1N N P DT W (d i , n). 1N N DT W (d i ) isthen one of the candidate considered by the heurisitic F DT W . Since F DT W returns the minimal classification error from all candidates, the classification error of 1N N DT W is always greater than or equal to F DT W . Algorithm 1 FDTW(training set, test set, n, nb rep=log(n)) # Look for a good value of the number of segments N # using the training set for (i in 0 : (nb rep -1)) do tab N ← 1 : (n -i) l ← f loor(n/sqrt(2 * n)) tab N candidats ← seq(f rom = n, to = 1, by = -l) # Parallel execution of 1NNPDTW mat r ← 1N N P DT W (training set, tab N candidats) # Mark candidates already used to not reuse for (i in tab candidats) do tab N [i] ← -1 end for # Search for the best candidate with the minimal error min ← minimun(mat r) # look for the local minimun near of the best candidate result[[(i + 1)]] ← localM inimun(min.N min, min.error min, training set, tab N ) end for # The best local minimal error m ← minimun(result) return m

Figure 1 :

 1 Figure 1: Comparison of the number of tested values of the parameter, number of segments with the Brute-force search algorithm, FDTW and IDDTW. x-axis datasets are sorted according to the length of the time series.

Figure 2 :

 2 Figure 2: Comparison of the execution time of the Brute-force search algorithm, FDTW and IDDTW.

Figure 3 :Figure 4 :

 34 Figure 3: Comparison of the classification error of the Brute-force search algorithm in x-axis and FDTW in y-axis. 1NN associated with Euclidiean distance, DTW without warping windows, DTW with warping windows. FDTW outperforms 1NN associated with Euclidiean distance

Table 1 :

 1 [START_REF] Bagnall | Time Se-ries Classification Website[END_REF]. In this paper we report the 1NN classification error of Bruteforce search, IDDTW and FDTW. Detailed information about the datasets 4.3 Discussions. Comparing FDTW and BF approaches (Figure

	N • Dataset	clas Size	Size	Time
			ses	of	of	se
				trai	tes	ries
				ning ting	len
				set	set	gth
	1	50Words	50	450	455	270
	2	Adiac	37	390	391	176
	3	Beef	5	30	30	470
	4	Car	4	60	60	577
	5	CBF	3	30	900	128
	6	Coffee	2	28	28	286
	7	Cricket X	12	390	390	300
	8	Cricket Y	12	390	390	300
	9	Cricket Z	12	390	390	300
	10	Distal	3	139	400	80
		Phalanx
		OutlineAge
		Group
	11	Distal	6	139	400	80
		Phalanx
		TW
	12	Earthquakes	2	139	322	512
	13	ECG	2	100	100	96
	14	ECGFiveDays 2	23	861	136
	15	Face (all)	14	560	1690	131
	16	Face (four)	4	24	88	350

Table 2 :

 2 Comparison of classification errors. In italics, the smallest classification error. In bold, the smallest classification error between IDDTW and FDTW. N is the number of segments selected and is the number of data points in a segment which is equal to n N .

	N • 1NN	1NN	1-NN	Smallest Brute N( )	IDDTW N( )	FDTW N( )
		Eucli	DTW DTW (r) known	force
		dean			error	search
		distance
	1	0.369	0.310	0.242 (6)	0.179	0.262	251(1)	0.268	256(1) 0.268	258(1)
	2	0.389	0.396	0.391 (3)	0.190	0.379	162(1)	0.432	128(1) 0.414	143(1)
	3	0.333	0.367	0.333 (0)	0.181	0.233	286(2)	0.3	8(59)	0.367	94(5)
	4	0.267	0.267	0.233 (1)	0.098	0.183	52(11)	0.367	8(72)	0.367	377(1)
	5	0.148	0.003	0.004 (11) 0.002	0	118(1)	0.003	128(1) 0.001	128(1)
	6	0.000	0.000	0.000 (0)	0	0	13(22)	0	64(4)	0.000	286(1)
	7	0.423	0.246	0.228 (10) 0.186	0.228	142(2)	0.256	256(1) 0.269	84(4)
	8	0.433	0.256	0.238 (17) 0.185	0.231	271(1)	0.241	256(1) 0.244	294(1)
	9	0.413	0.246	0.254 (5)	0.173	0.221	249(1)	0.223	256(1) 0.233	276(1)
	10	0.218	0.208	0.228 (1)	0.171	0.2	78(1)	0.225	16(5)	0.223	80(1)
	11	0.273	0.29	0.272 (0)	0.272	0.263	35(2)	0.288	16(5)	0.278	80(1)
	12	0.326	0.258	0.258 (22) 0.241	0.198	176(2)	0.258	512(1) 0.276	101(5)
	13	0.120	0.230	0.120 (0)	0.110	0.13	38(3)	0.19	8(12)	0.180	11(9)
	14	0.203	0.232	0.203 (0)	0.014	0.117	11(12)	0.289	32(4)	0.117	11(12)
	15	0.286	0.192	0.192 (3)	0.010	0.091	79(2)	0.194	128(1) 0.148	99(1)
	16	0.216	0.170	0.114 (2)	0.004	0.08	107(3)	0.352	32(11) 0.102	140(3)
	17	0.217	0.177	0.154(4)	0.026	0.154	149(3)	0.257	16(29) 0.177	27(17)
	18	0.087	0.093	0.087 (0)	0.001	0.02	38(4)	0.073	32(5)	0.020	38(4)
	19	0.4	0.533	0.400 (0)	0.164	0.343	21(20)	0.026	32(13) 0.432	32(13)
	20	0.630	0.623	0.588 (2)	0.483	0.549	328(3)	0.588	64(17) 0.594	948(1)
	21	0.658	0.616	0.613 (14) 0.448	0.578	1770(1) 0.627	64(29) 0.622	171(11)
	22	0.045	0.050	0.045 (0)	0.030	0.033	20(1)	0.043	8(3)	0.033	24(1)
	23	0.246	0.131	0.131 (6)	0.131	0.082	70(9)	0.246	16(40) 0.180	524(1)
	24	0.425	0.274	0.288 (5)	0.201	0.192	150(2)	0.301	64(5)	0.301	51(6)
	25	0.316	0.263	0.253 (20) 0.215	0.255	95(1)	0.271	64(2)	0.280	34(3)
	26	0.26	0.25	0.253 (5)	0.178	0.233	27(2)	0.283	2(40)	0.283	2(40)
	27	0.439	0.416	0.419 (2)	0.413	0.398	27(2)	0.414	32(2)	0.416	80(1)
	28	0.121	0.165	0.134 (1)	0.083	0.135	14(6)	0.197	16(5)	0.165	84(31)
	29	0.133	0.167	0.133 (0)	0.099	0.1	191(3)	0.167	32(18) 0.100	234(2)
	30	0.479	0.409	0.388 (7)	0.033	0.364	31(14)	0.372	32(13) 0.409	35(12)
	31	0.038	0	0.000 (6)	0	0	35(4)	0	128(1) 0	135(1)
	32	0.292	0.263	0.263 (6)	0.185	0.24	75(1)	0.288	4(20)	0.288	4(20)
	33	0.461	0.35	0.300 (3)	0	0.083	54(9)	0.239	64(7)	0.122	48(10)
	34	0.305	0.275	0.305 (0)	0.082	0.206	37(2)	0.208	16(4)	0.304	26(3)
	35	0.141	0.169	0.141 (0)	0.04	0.14	5(13)	0.197	16(4)	0.178	45(1)
	36	0.211	0.208	0.154 (2)	0.033	0.165	59(2)	0.195	64(2)	0.208	55(2)
	37	0.100	0.050	0.062 (8)	0.037	0.044	376(1)	0.059	32(12) 0.060	34(12)
	38	0.120	0.007	0.017 (6)	0.001	0.007	60(1)	0.437	2(30)	0.007	60(1)
	39	0.240	0.000	0.010 (3)	0	0	47(6)	0	64(4)	0.000	275(1)
	40	0.090	0.000	0.002 (4)	0	0	21(6)	0	64(2)	0.000	128(1)
	41	0.253	0.096	0.132 (5)	0.006	0.045	55(1)	0.073	32(3)	0.112	70(1)
	42	0.005	0.020	0.005 (1)	0	0.007	109(1)	0.013	8(19)	0.008	95(2)
	43	0.389	0.426	0.389 (0)	0.074	0.204	3(78)	0.463	20(11) 0.37	128(1)
	44	0.382	0.351	0.252 (8)	0.222	0.337	133(2)	0.365	64(4)	0.343	135(2)
	45	0.170	0.164	0.155 (2)	0.090	0.149	117(4)	0.158	128(3) 0.154	384(1)
	X	0.268	0.227	0.242		0.175		0.214		0.232


Ariela Briani

Pierre Cardaliaguet
email: cardaliaguet@ceremade.dauphine.fr

Stable solutions in potential mean field game systems

We introduce the notion of stable solution in mean field game theory: they are locally isolated solutions of the mean field game system. We prove that such solutions exist in potential mean field games and are local attractors for learning procedures.

Introduction

Mean field games (MFG) are Nash equilibrium configurations in differential games with infinitely many infinitesimal players. If the existence of such equilibria holds in general frameworks, one cannot expect their uniqueness without strong restrictions. However the multiplicity of equilibria is a real issue in terms of applications: indeed, when there are several equilibria, it is not clear how the players can coordinate in order to decide which equilibrium to play. The aim of this work is to introduce a notion of stable solutions for the MFG system: these stable solutions are Nash equilibria which have the property of being locally isolated as well as local attractors for learning procedures. In this sense they are robust with respect to perturbation.

Let us recall that the terminology and main properties of Mean Field Games were introduced by Lasry and Lions in a series of papers [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]. At the same period Huang, Caines and Malhamé [START_REF] Huang | Large population stochastic dynamic games: closed-loop McKean-Vlasov systems and the Nash certainty equivalence principle[END_REF] discussed the same concept under the name of Nash certainty equivalence principle. One way to represent an MFG equilibrium is through the following system of partial differential equations.

   -∂ t u -∆u + H(x, Du) = f (x, m) in T d × (t 0 , T ), ∂ t m -∆m -div(mD p H(x, Du)) = 0 in T d × (t 0 , T ) m(x, t 0 ) = m 0 (x), u(x, T ) = g(x, m(T )) in T d .

(1) (To simplify the discussion related to boundary conditions, we work here in the torus T d := R d /Z d ). In the above system the pair (u, m) is the unknown. The map u = u(x, t) can be interpreted as the value function of a (small) player, while m(t) = m(x, t) is understood as the evolving probability density of the players at time t. Note that u satisfies a backward Hamilton-Jacobi equation while m solves a forward Kolmogorov equation with initial condition m(t 0 ) = m 0 , where m 0 is a given probability density on T d . The Hamiltonian H : T d × R d → R is typically smooth and convex in the second variable. The coupling functions f, g depend on the space variable and on the probability density.

Existence of solution for the MFG system (1) has been established under general assumptions by Lasry and Lions [START_REF] Lasry | Jeux à champ moyen. II. Horizon fini et contrôle optimal[END_REF][START_REF] Lasry | Mean field games[END_REF]. Uniqueness, however, only holds under restrictive conditions: either the horizon is small, or the couplings functions are monotone, see [START_REF] Lasry | Mean field games[END_REF].

In this paper we are interested in MFG systems which might have several solutions. Our aim is to introduce a particular notion, the stable MFG equilibrium. The natural idea for this is to call stable a solution of the MFG system [START_REF] Ambrosio | Gradient flows in metric spaces and in the space of probability measures[END_REF] for which the associated linearized system has only one solution, the trivial one. We first show (Proposition 4.2) that, with this definition, stable MFG solutions are isolated, in the sense that there is no other solution (with the same initial condition) in a neighborhood. If this result is very natural and quite expected, it illustrates well the notion.

The two other results are more subtle and obtained in the framework of potential MFG games. Following [START_REF] Cardaliaguet | Learning in Mean Field Games: the Fictitious Play[END_REF], we say that the MFG game is potential if there exists F, G such that δF δm (m, x) = f (x, m) and δG δm (m, x) = g(x, m),

where the derivative (with respect to the measure m) in taken in the sense of [START_REF] Cardaliaguet | The master equation and the convergence problem in mean field games[END_REF] (see also subsection 2.1). Following [START_REF] Lasry | Mean field games[END_REF], we know then that the MFG game has a "potential". Namely, let us define the functional (where m = m(x, t) is an evolving probability density and w = w(x, t) is a vector field)

J(m, w) := T t0 T d L x, w (x, t) m(x, t) m(x, t)dxdt +
T t0 F (m(t))dt + G(m(T ))
where L is the convex conjugate of H (see [START_REF] Cannarsa | Semiconcave functions, Hamilton-Jacobi equations and optimal control[END_REF]) and the pair (m, w) has to fulfill the continuity equation

∂ t m -∆m + div(w) = 0 in T d × [t 0 , T ], m(t 0 ) = m 0 .
Then any minimizer (m, w) of J corresponds to a solution to the MFG system [START_REF] Lasry | Mean field games[END_REF], in the sense that there exists u = u(x, t) such that the pair (u, m) solves (1) and w = -mD p H(x, Du).

Our first main result says that, when the game is potential, there are many initial measures m 0 starting from which there is a stable solution to the MFG system. Indeed, we show (Theorem 4.3) that, if (u, m) is a MFG equilibrium on the time interval [t 0 , T ] which corresponds to a minimizer of the potential J, then its restriction to any subinterval [t 1 , T ] (where t 1 ∈ (t 0 , T )) is a stable MFG equilibrium.

Second we show (in Theorem 5.1) that, for potential MFG systems, stable equilibria are local attractors for a learning procedure. Here we consider the learning procedure inspired by the Fictitious Play [START_REF] Brown | Iterative solution of games by Fictitious Play[END_REF] and introduced by the second author and S. Hadikhanloo in [START_REF] Cardaliaguet | Learning in Mean Field Games: the Fictitious Play[END_REF]. Given µ 0 = (µ 0 (t)) an initial guess of the evolving probability density of the players, we define by induction the sequence (u n , m n , µ n ):

   -∂ t u n+1 -∆u n+1 + H(x, Du n+1 ) = f (x, µ n ) ∂ t m n+1 -∆m n+1 -div(m n+1 D p H(x, Du n+1 )) = 0 m n+1 (0) = m 0 , u n+1 (T, x) = g(x, µ n (T ))
and

µ n+1 = n n + 1 µ n + 1 n + 1 m n+1 . ( 2
)
The interpretation of this system is that the game is played over and over. At stage n + 1 (the (n + 1)-th time the game is played), all the players play as if the population density is going to evolve according to µ n . The map u n+1 is the value function of each small player. When the players play optimally in this game, the population density actually evolves according to m n+1 . The players then update their estimate of the evolving population density by taking the average of the previous observations (i.e., by the rule (2)). It is proved in [START_REF] Cardaliaguet | Learning in Mean Field Games: the Fictitious Play[END_REF] that this simple learning procedure converges: namely, if the game is potential, then any converging subsequence of the relatively compact (for the uniform convergence) sequence (u n , m n ) is a MFG Nash equilibrium.

Here we show that, if in addition the equilibrium (u, m) is stable and if µ 0 is sufficiently close to m, then the entire sequence (u n , m n ) converges to (u, m). This result illustrates again the robustness of the stable equilibria: if the players deviate from a stable equilibrium configuration, the learning procedure pushes them back to this equilibrium.

The techniques used in the paper are inspired by finite dimensional optimal control (see, for instance, in the monograph of Cannarsa and Sinestrari [START_REF] Cannarsa | Semiconcave functions, Hamilton-Jacobi equations and optimal control[END_REF]). In particular the fact that the solutions are stable when restricted to a subinterval is known in this context. The proof requires a uniqueness result for the solution of the MFG system and its linearized version given an initial condition for m and for Du. To show such a statement we rely on a method developed by Lions and Malgrange for the backward uniqueness of the heat equation, and subsequently extended to systems in Cannarsa and Tessitore [START_REF] Cannarsa | Optimality conditions for boundary control problems of parabolic type[END_REF].

The paper is organized in the following way: we first introduce the notation, assumptions and recall standard existence and uniqueness results for the MFG system. We also prove a first new uniqueness result for the MFG system given the initial measure and the initial vector field. Then we discuss the notion of potential games (section 3) and illustrate the notion by expliciting an example of multiple solutions for a MFG system. In section 4 we define the notion of stable MFG equilibria and provide our first main result (Theorem 4.3) on the existence of such equilibria in the potential case. We complete the paper by the analysis of the fictitious play for MFG system (section 5). In Appendix we prove the uniqueness of a general linear forward-backward system with given initial data: this result is used several times in the text.

Acknowledgement:

The authors were partially supported by the ANR (Agence Nationale de la Recherche) project ANR-16-CE40-0015-01.

2 Assumptions and basic results on MFG systems

Notation and assumption

Throughout the paper we work on the d-dimensional torus T d := R d /Z d : this simplifying assumption allows us to ignore issues related to boundary conditions. We also work in a finite horizon T > 0. We denote by P(T d ) the set of Borel probability measures on T d , endowed with the Monge-Kantororitch distance d 1 . We define M(T d , R d ) as the set of Borel vector measures w with finite mass |w|. For α ∈ [0, 1], we denote by C α ([0, T ], P(T d )) the set of maps m : [0, T ] → P(T d ) which are α-Holder continuous if α ∈ (0, 1), continuous if α = 0, Lipschitz continuous if α = 1.

Next we recall the notion of derivative of a map U : P(T d ) → R as introduced in [START_REF] Cardaliaguet | The master equation and the convergence problem in mean field games[END_REF]. We say that U is C 1 if there exists a continuous map δU δm :

T d × P(T d ) → R such that U (m ′ ) -U (m) = 1 0 T d δU δm (x, (1 -t)m + tm ′ )(m ′ -m)(dx)dt ∀m, m ′ ∈ P(T d ).
The derivative is normalized by the condition

T d δU δm (x, m)dm(x) = 0 ∀m ∈ P(T d ).
We write indifferently δU δm (x, m)(µ) and

T d δU δm (x, m)dµ(x)
for a signed measure µ with finite mass.

If u : T d × [0, T ] → R is a sufficiently smooth map, we denote by Du(x, t) and ∆u(x, t) its spatial gradient and spatial Laplacian and by ∂ t u(x, t) its partial derivative with respect to the time variable. For p = 1, 2, ∞, we denote by • p the L p norm of a map on T d . We denote by C 0 the set of continuous maps, by C 2,0 the set of of maps such that D 2 u and ∂ t u exist and are continuous. By abuse of notation, we set

u C 1,0 = u ∞ + Du ∞ .
We will also use the classical Holder space. For α ∈ (0, 1), we denote by C 0,α the set of map u = u ∈ C 0 which are α-Holder continuous in space and α/2 in time. The set C 1,α is the set of maps u ∈ C 0 such that u and Du belong to C 0,α . Finally C 2,α consists in the maps u ∈ C 2,0 such that D 2 u and ∂ t u belong to C 0,α . Let us recall that, if u is in C 2,α , then u is also in C 1,α .

Assumptions. The following assumptions are in force throughout the paper.

• The Hamiltonian H = H(x, p) : T d × R d → R is of class C 2 and satisfies C -1 I d ≤ D 2 pp H(x, p) ≤ CI d (3)
We define the convex conjugate L of H as L(x, q) = sup

p∈R d {-p • q -H(x, p)}. (4)
• The coupling functions f, g :

T d × P(T d ) → R

The MFG system

Let us recall that, under our standing assumptions, the MFG system (1) has at least one classical solution: see, for instance, [START_REF] Lasry | Mean field games[END_REF]. In general it is not expected that this solution is unique: uniqueness in known to hold for short time horizon (or small data) or when the coupling functions f and g are monotone [START_REF] Lasry | Mean field games[END_REF]. We provide in the next section an example of multiple solutions. However we prove here that there is only one solution given the initial measure and the initial vector field. This result, which is of limited interest for the true mean field game system, will be used several times in the sequel. Proposition 2.1. Let (u 1 , m 1 ) and (u 2 , m 2 ) be two solutions of the MFG system (1) with the same initial initial condition for the measure m 1 (t 0 ) = m 2 (t 0 ) = m 0 and the same initial vector field Du

1 (•, t 0 ) = Du 2 (•, t 0 ) in T d . Then (u 1 , m 1 ) = (u 2 , m 2 ) on T d × [t 0 , T ].
Proof. Without loss of generality we assume that t 0 = 0. Let us set, for i = 1, . . . , d,

w i = ∂ xi (u 2 -u 1 ), w = (w i ) i=1,...,d , µ = m 2 -m 1 . Then (w, µ) solves the system    -∂ t w i -∆w i + g i = 0 in T d × (0, T ), i = 1, . . . , d, ∂ t µ -∆µ + div(h) = 0 in T d × (0, T ) w i (x, 0) = µ(x, 0) = 0 in T d , i = 1, . . . , d.
where

g i (x, t) = ∂ xi H(x, Du 2 ) + D p H(x, Du 2 ) • D(∂ xi u 2 ) -∂ xi H(x, Du 1 ) -D p H(x, Du 1 ) • D(∂ xi u 1 ) + ∂ xi f (x, m 2 ) -∂ xi f (x, m 1 ) and h(x, t) = m 2 D p H(x, Du 2 (x, t)) -m 1 D p H(x, Du 1 (x, t)) .
As the u i are classical solutions and, by our assumption, ∂ xi f :

T d × P(T d ) → R is globally
Lipschitz continuous, we have

|g i (x, t)| ≤ C (|D(u 2 -u 1 )(x, t)| + |D(∂ xi (u 2 -u 1 ))(x, t)| + m 2 (•, t) -m 1 (•, t) 2 ) ≤ C (|w(x, t)| + |Dw(x, t)| + µ(•, t) 2 ) .
In the same way,

|h(x, t)| ≤ C (|(m 2 -m 1 )(x, t)| + |D(u 2 -u 1 )(x, t)|) ≤ C (|w(x, t)| + |µ(x, t)|) . Moreover, as Du i ∞ + D 2 u i ∞ + m i ∞ + Dm i ∞ ≤ C, we also obtain |div(h)(x, t)| ≤ C (|w(x, t)| + |Dw(x, t)| + |µ(x, t)| + |Dµ(x, t)|) .
Then we can conclude from Theorem 6.1 in the Appendix that µ = 0 on T d × [0, T ]. Since m 1 = m 2 and since the backward Hamilton-Jacobi equations for the u i have the same right-hand side, their solution is the same and we conclude that u 1 = u 2 .

Potential MFG systems

Throughout this section we assume that the MFG system is potential, i.e., the coupling functions f and g derive from potentials: there exists F, G :

P(T d ) → R of class C 1 such that f = δF δm , g = δG δm . (5)
The aim of the section is to recall that, under this condition, the MFG systems is then potential (i.e., the equilibria are obtained through a minimization procedure). We derive from this an example of non uniqueness of the solution of the MFG system.

Solutions of the MFG system and minimizers of the potential

Let t 0 ∈ [0, T ] and m 0 ∈ P(T d ). Our aim is to define the functional

J J(m, w) = T t0 T d L x, w(x, t) m(x, t) m(x, t)dxdt + T t0 F (m(t))dt + G(m(T ))
under the constraint

∂ t m -∆m + div(w) = 0 in T d × [t 0 , T ], m(t 0 ) = m 0 ,
for a large class of data (m, w). This classical construction is reminiscent of the Benamou and Brenier approach of optimal transport and can be found, for instance, in [START_REF] Ambrosio | Gradient flows in metric spaces and in the space of probability measures[END_REF]. We denote by E(t 0 ) the set of pairs time dependent Borel measures (m(t), w(t)

) ∈ P(T d ) × M(T d , R d ) such that t → m(t) is continuous, T 0 |w(t)|dt < ∞,
and equation

∂ t m -∆m + div(w) = 0 in T d × [t 0 , T ], m(t 0 ) = m 0
holds in the sense of distribution. We also denote by E 2 (t 0 ) the subset of (m(t), w(t)) ∈ E(t 0 ) such that w(t) is absolutely continuous with respect to m(t) with a density dw(t) dm(t) satisfying

T d T 0 dw(t) dm(t) (x) 2 m(dx, t)dt < ∞.
Then we defined J on E(t 0 ) by

J(m, w) :=            T t0 T d L x, dw(t) dm(t) (x) m(dx, t)dt + T t0 F (m(t))dt + G(m(T )) if (m, w) ∈ E 2 (t 0 )
+∞ otherwise When it will be important to stress the dependence on the initial data we will write J(t 0 , m 0 , •, •).

The following result states that minimizers of the functional J correspond to solution of the MFG system. This remark was first pointed out in [START_REF] Lasry | Mean field games[END_REF] and used repetitively since then in different contexts.

Proposition 3.1. Under our standing assumptions:

(i) For any t 0 ∈ [0, T ] and m 0 ∈ P(T d ) there exists a minimum (m, w) ∈ E 2 of J(t 0 , m 0 , •, •).

(ii) Let (m, w) be minimum of J(t 0 , m 0 , •, •). Then there exists u such that (u, m) is a classical solution to the MFG system

   -∂ t u -∆u + H(x, Du) = f (x, m) in T d × (t 0 , T ), ∂ t m -∆m -div(mD p H(x, Du)) = 0 in T d × (t 0 , T ) m(x, t 0 ) = m 0 (x), u(x, T ) = g(x, m(T )) in T d (6)
and w(x, t) = -m(x, t)D p H(x, Du). In particular, any minimizer is a classical solution of the above system, i.e., u, m ∈ C 2,α .

The proof is standard and is closely related to techniques used in optimal transport theory. However, as it has never been explicitly checked in this specific framework where the potentials are non convex, we provide the main argument for the sake of completeness. Proof. To prove (i) let us consider a minimizing sequence (m n , w n ) ∈ E 2 . By construction J(m n , w n ) ≤ C, thus the coercivity assumption (3) on H implies the following uniform bound:

T t0 T d dw n (t) dm n (t) (x) 2 m n (dx, t)dt ≤ C . (7)
We can then argue as in [START_REF] Cardaliaguet | Second order mean field games with degenerate diffusion and local coupling[END_REF]Lemma 3.1] to conclude that the sequence (m n ) is uniformly bounded in C 1/2 ([0, T ], P(T d )). In particular,

T t0 |w n (t)|dt ≤ T t0 T d dw n (t) dm n (t) 2 m n (dx, t)dt 1/2 T t0 T d m n (dx, t)dt 1/2 ≤ C.
So, up to a subsequence,

w n → w in M((0, T ) × T d , R d ) and (m n ) converges in C 0 ([0, T ], P(T d )).
By standard argument the pair (m, w) belongs to E 2 . To conclude that the couple (m, w) is a minimum of the functional J we remark that the functional J is lower semicontinuous on E.

In order to prove (ii) we define on E

Φ(m, w) :=        T t0 T d L x, w(t) m(t) m(dx, t)dt if (m, w) ∈ E 2 +∞ otherwise and Ψ(m) := T t0 F (m(t))dt + G(m(T )). Note that J(m, w) = Φ(m, w) + Ψ(m). Let ( m, w) be a minimum of J. Recall that m ∈ C 1/2 ([0, T ], P(T d ))
. We first claim that, for any (m, w) ∈ E, we have

- T t0 T d f (x, m(t))(m(dx, t) -m(dx, t)) - T d g(x, m(T ))(m(dx, T ) -m(dx, T )) ≤ Φ(m, w) -Φ( m, w) . (8)
Indeed, setting m λ := (1 -λ) m + λm, w λ := (1 -λ) w + λw , λ ∈ (0, 1) we have by minimality of ( m, w):

Φ(m λ , w λ ) -Φ( m, w) ≥ Ψ( m) -Ψ(m λ ).
Thus, by the regularity assumptions on F and G and the convexity of Φ, we obtain

λ(Φ(m, w) -Φ( m, w)) ≥ λ - T t0 T d f (x, m(t, x))(m -m)(dx, t) - T d g(x, m(x, T ))(m -m)(dx, T )) + o(λ).
Then [START_REF] Cardaliaguet | The master equation and the convergence problem in mean field games[END_REF] follows dividing by λ and letting λ tends to 0. By [START_REF] Cardaliaguet | The master equation and the convergence problem in mean field games[END_REF], the pair ( m, w) is a minimizer of the following (local and convex) functional on E:

J(m, w) = Φ(m, w) + T 0 T d f (x, m(t))m(dx, t)dt + T d g(x, m(T ))m(dx, T ).
We can then follow the standard arguments for convex functionals (see for instance [START_REF] Cardaliaguet | Second order mean field games with degenerate diffusion and local coupling[END_REF]): the problem of minimizing J on E is the dual problem (in the sense of the Fenchel-Rockafellar duality theorem) of the problem

inf u∈C 2 - T d m 0 (x)u(x, 0)dx : -∂ t u-∆u+H(x, Du) ≤ f (x, m(t)) and u(x, T ) ≤ g(x, m(x, T )) .
By comparison, there is an obvious minimum to this problem which is the solution ū to

-∂ t ū -∆ū + H(x, D ū) = f (x, m(t)) in T d × (0, T ) u(x, T ) = g(x, m(x, T )) in T d . This solution is C 2,α because m ∈ C 1/2 ([0, T ), P(T d ))
. By the Fenchel-Rockafellar duality theorem, see also [6, Lemma 3.2], we have that 0 = J( m, w) -

T d m 0 (x)ū(x, 0)dx .
This means that

0 = T t0 T d L x, dw(t) dm(t) + f (x, m(t)) m(dx, t)dt + T d g(x, m(T )) m(dx, T ) - T d m 0 (x)ū(x, 0)dx = T t0 T d L x, dw(t) dm(t) + (-∂ t ū -∆ū + H(x, D ū)) m(dx, t)dt + T d g(x, m(T )) m(dx, T ) - T d m 0 (x)ū(x, 0)dx = T t0 T d L x, dw(t) dm(t) + H(x, D ū) + Du • ( d w(t) d m(t) ) m(dx, t)dt
where we used the equation for ū in the second equality and the equation for ( m, w) in the last one. Recalling that L is the convex conjugate of H which is uniformly convex, we find

d w(t) d m(t) = -D p H(x, D ū) m -a.e.
This means that m solves the Kolmogorov equation

∂ t m -∆ m -div( mD p H(x, D ū)) = 0, m(•, t 0 ) = m 0 ,
which has a regular drift: thus m is of class C 2,α by Schauder theory. Therefore w is also smooth and the proof of (ii) is complete.

We now explain that, if (m, w) is a minimizers of J(t 0 , m 0 , •, •), then, for any later time t 1 > t 0 , (m, w) is the unique minimizer for J(t 1 , m(t 1 ), •, •). The argument borrows ideas in finite dimensional control theory [START_REF] Cannarsa | Semiconcave functions, Hamilton-Jacobi equations and optimal control[END_REF]. Proposition 3.2. Let (m, w) be minimum of J(t 0 , m 0 , •, •) for a given initial condition m 0 at time t 0 . Then, for any t 1 ∈ (t 0 , T ), J(t 1 , m(t 1 ), •, •) has a unique minimum, which is the restriction to [t 1 , T ] of (m, w).

Proof. By dynamic programming principle, the restriction to

[t 1 , T ] of (m, w) is a minimum of J(t 1 , m(t 1 ), •, •). Let ( m, w) be another minimum of J(t 1 , m(t 1 ), •, •). Then the following map ( m, ŵ) = (m, w) on T d × [t 0 , t 1 ) ( m, w) on T d × [t 1 , T ]
is also a minimum of J(t 0 , m 0 , •, •). Thus, by Proposition 3.1, there exist u, û such that (u, m) and (û, m) are both classical solutions to the MFG system [START_REF] Cardaliaguet | Second order mean field games with degenerate diffusion and local coupling[END_REF], with w = -mD p H(•, Du) and ŵ = -mD p H(•, D û). On T d × (t 0 , t 1 ), we have m = m > 0 and -mD p H(•, D û) = -mD p H(•, Du), so that, by uniform convexity of H, we obtain the equality Du = D û on T d × (t 0 , t 1 ). By continuity we infer that Du(•, t

1 ) = D û(•, t 1 ). Then Proposition 2.1 implies that (û, m) = (ũ, m) = (u, m) on T d × [t 1 , T ].

An example of multiple solutions

As an application of Proposition 3.1, we provide here a simple example of multiple solutions for the MFG system. Such examples in the time dependent setting are scarce: we are only aware of a recent result by Bardi and Fischer [2]. The idea is relatively elementary: we suppose that the MFG system is symmetric with respect to the x 1 variable, but that the potential favors asymmetric solutions. Then the MFG system enjoys a symmetric solution as well as an asymmetric one.

To build the example, let τ (x 1 , x 2 , . . . , x d ) = (-x 1 , x 2 , . . . , x d ) (viewed as a map on R d as well as on T d ) and let us denote by P s (T d ) the subset of measures m ∈ P(T d ) such that τ ♯m = m. We suppose, in addition to our general conditions on H, f and g, that:

(1) The function F is nonnegative and symmetric with respect to the x 1 variable:

F (τ ♯ m) = F (m), ∀m ∈ P(T d ).
Note that this implies that f (τ (x), m) = f (x, m) for any m ∈ P s (T d ).

(2) There exists a measure µ ∈ P(T d ) such that inf

m∈Ps(T d ) F (m) > F (µ).
(3) The Hamiltonian H is symmetric with respect to the x 1 variable:

H(τ (x), τ (p)) = H(x, p) for any (x, p) ∈ T d × R d . Moreover L is nonnegative.
Given a parameter θ > 0, we consider the MFG system

   -∂ t u -∆u + H(x, Du) = θ f (x, m) in T d × (0, T ), ∂ t m -∆m -div(mD p H(x, Du)) = 0 in T d × (0, T ) m(x, 0) = m 0 (x), u(x, T ) = 0 in T d . (9)
Proposition 3.3. Under our assumptions, for any symmetric initial condition m 0 ∈ P s (T d ) with a smooth density, there exist T > 0 and θ > 0 large enough such that the MFG system (9) has at least two different solutions.

Proof. Let us first explain that the MFG system has a symmetric solution. For this we build, in a standard way, a fixed point mapping Φ on C α ([0, T ], P s (T d )) (for some fixed α ∈ (0, 1/2)): see for instance [START_REF] Lasry | Mean field games[END_REF]. To m ∈ C α ([0, T ], P s (T d )) we first associate the solution to

-∂ t u -∆u + H(x, Du) = θ f (x, m) in T d × (0, T ), u(x, T ) = 0 in T d .
Because of our assumptions, the above equation has a unique classical solution u, which is therefore symmetric with respect to x 1 . Moreover, thanks to the uniform convexity of H, Du is bounded independently of m. Next we consider the solution m to the Kolmogorov equation

∂ t m -∆ m -div( mD p H(x, Du)) = 0 in T d × (0, T ) m(x, 0) = m 0 (x) in T d .
Here again there exists a classical solution m, which is symmetric with respect to x 1 . Moreover, as Du is bounded, m is bounded in C 1/2 ([0, T ], P s (T d )) independently of m. Finally we set Φ(m) := m. It is easy to check that Φ has a fixed point, which leads to a symmetric solution to the MFG system. We denote by (û, m) this symmetric solution and set ŵ := -mD p H(x, D û).

Next we show that the symmetric solution cannot minimize the functional J for suitable choice of T and θ. Let us denote by M the minimum value of the functional F over P s (T d ).

Recalling that there exists µ such that F (µ) < M . We fix ε, η > 0 such that

M > ε + max d1(m,µ)≤η F (m).
We then construct (e.g. by convolution) a measure m ∈ P(T d ) which is absolutely continuous and such that (for some constant

C > 0) m ∈ C ∞ (T d ) , d 1 ( m, µ) ≤ η , 1 C ≤ m(x) ≤ C.
Setting w = D m, we have -∆ m -div( w) = 0. Then we choose θ ≥ 1 large enough to have

εθ 2 > T d L x, w(x) m(x) m(x)dx.
Next we define the pair (m, w) which connects m 0 and m in time 1 and is equal to ( m, w) after time 1:

(m(x, t), w(x, t)) := ((1 -t)m 0 (x) + t m(x), -Dφ(x) + (1 -t)Dm 0 (x) + tD m(x)) for t ∈ [0, 1],
(where φ is a solution to -∆φ = m 0 -m in T d ) and

(m(x, t), w(x, t)) = ( m(x), w(x)) for t ∈ [1, T ].
Then one easily checks that the pair (m, w) belongs to E 2 . Finally we choose T large enough to have sup

m ′ ∈P(T d ) F (m ′ ) + 1 0 T d m(x, t)L x, w(x, t) m(x, t) dxdt ≤ εθT 2 .
As L ≥ 0, we have

J( m, ŵ) ≥ θ T M > θT ε + θT max d1(m ′ ,µ)≤η F (m ′ ).
On the other hand,

J(m, w) = 1 0 T d L(x, w m )mdx + F (m(t)) dt + (T -1) T d L(x, w m ) mdx + F ( m) < εθT /2 + (T -1)εθ/2 + (T -1) max d1(m ′ ,µ)≤η F (m ′ ) < J( m, ŵ).
This proves that ( m, ŵ) is not a minimizer of J. As, by Proposition 3.1, J has a minimum which is associated with a solution of the MFG system (9), there exists a solution to the MFG system (9) different from (û, m).

Stable solution to the MFG system 4.1 Definition and basic property

We say that a solution of the MFG system (1) is stable if the unique solution to linearized system is the trivial one. Definition 4.1. Let (u, m) be a solution to the MFG system (1) with initial condition (t 0 , m 0 ) ∈ [0, T ] × P(T d ). We say that the solution (u, m) is stable if (v, µ) = (0, 0) is the unique solution to the linearized system

   -∂ t v -∆v + D p H(x, Du) • Dv = δf δm (x, m)(µ) in T d × (t 0 , T ), ∂ t µ -∆µ -div(µD p H(x, Du)) -div(mD 2 pp H(x, Du)Dv) = 0 in T d × (t 0 , T ) µ(x, t 0 ) = 0, v(x, T ) = δg δm (x, m(T ))(µ(T )) in T d .
Before proving that stable solutions do exist, we show that they are isolated.

Proposition 4.2. Let (u, m) be a stable solution starting from an initial position (t 0 , m 0 ) ∈ [0, T ] × P(T d ). Then there is η > 0 such that, for any m 1 0 ∈ P(T d ) with m 1 0 -m(t 0 ) C 0 ≤ η, there exists at most one solution (u 1 , m 1 ) of the MFG system

   -∂ t u 1 -∆u 1 + H(x, Du 1 ) = f (x, m 1 ) in T d × (t 0 , T ), ∂ t m 1 -∆m 1 -div(m 1 D p H(x, Du 1 )) = 0 in T d × (t 0 , T ) m 1 (x, t 0 ) = m 1 0 (x), u 1 (x, T ) = g(x, m 1 (T )) in T d . such that (u, m) -(u 1 , m 1 ) C 1,0 ×C 0 ≤ η.
Proof. We argue by contradiction, assuming that there exists two distinct solutions (u n,1 , m n,1 ) and (u n,2 , m n,2 ) of the MFG system with m n,1 (•, t 0 ) = m n,2 (•, t 0 ) and converging to (u, m) as

n → +∞ in C 1,0 × C 0 . We set (v n , µ n ) = (ρ n ) -1 (u n,2 , m n,2 ) -(u n,1 , m n,1 )
,

where ρ n := (u n,2 , m n,2 ) -(u n,1 , m n,1 ) C 1,0 ×C 0 .
We note that the pair (v n , µ n ) solves

-∂w -∆w + g n = 0 in T d × (t 0 , T ), ∂µ -∆µ + div(h n ) = 0 in T d × (t 0 , T ),
where

g n (x, t) = (ρ n ) -1 H(x, Du n,2 (x, t)) -H(x, Du n,1 (x, t)) -f (x, m n,2 (t)) + f (x, m n,1 (t))
and

h n (x, t) = (ρ n ) -1 -m 2,n (x, t)D p H(x, Du n,2 (x, t)) + m 1,n (x, t)D p H(x, Du n,1 (x, t)) .
Note that g n ∞ and h n ∞ are bounded by definition of ρ n . Thus, by standard arguments in uniform parabolic equations (see for instance, Theorem 9.1 and Theorem 10.1 in [START_REF] Ladyženskaja | N Linear and quasilinear equations of parabolic type[END_REF]), there exists α ∈ (0, 1) such that (v n ) is bounded in C 1+α,(1+α)/2 while (µ n ) is bounded in C α,α/ε . Thus, up to a subsequence denoted in the same way, (v n , µ n ) converges in C 1,0 × C 0 to some (v, µ). Note that (v, µ) C 1,0 ×C 0 = 1 and that (v, µ) is a (a priori weak) solution to the linearized problem

   -∂ t v -∆v + D p H(x, Du) • Dv = δf δm (x, m)(µ) in T d × (t 0 , T ), ∂ t µ -∆µ -div(µD p H(x, Du)) -div(mD 2 pp H(x, Du)Dv) = 0 in T d × (t 0 , T ) µ(x, t 1 ) = 0, v(x, T ) = δg δm (x, m(T ))(µ(T )) in T d .
By parabolic regularity, (v, µ) is actually a classical solution to the above equation. As (u, m) is a stable solution, one must have (v, µ) = (0, 0). This leads to a contradiction with the fact that (v, µ) C 1,0 ×C 0 = 1.

Existence of stable solutions for potential MFG systems

The next result states that, if (u, m) is a solution to a potential MFG system on the time interval [t 0 , T ] which corresponds to a minimizer of J, then the restriction of (u, m) to any subinterval [t 1 , T ] (where t 1 ∈ (t 0 , T )) is a stable solution of the MFG system.

Theorem 4.3. Let us assume that the game is potential, i.e, satisfies (5) for some C 1 maps F, G : P(T d ) → R. Let (m, w) be minimum of J(t 0 , m 0 , •, •) for a given initial condition m 0 ∈ P(T d ) at time t 0 ∈ [0, T ]. Let u be such that (u, m) is a solution to the MFG system (1). Then, for any t 1 ∈ (t 0 , T ), the restriction of (u, m) to the time interval [t 1 , T ] is a stable solution to the MFG system.

Note that the Theorem implies that there are many stable solutions to the MFG system. In particular, the set of initial conditions (t 0 , m 0 ) for which there is a stable solution is dense.

We now start the proof of Theorem 4.3. We have to show that, for any t 1 ∈ (t 0 , T ), (v, µ) := (0, 0) is the unique classical solution to the linearized system

         -∂ t v -∆v + D p H(x, Du) • Dv = δf δm (x, m)(µ) in T d × (t 1 , T ), ∂ t µ -∆µ -div(µD p H(x, Du)) -div(mD 2 pp H(x, Du)Dv) = 0 in T d × (t 1 , T ) µ(x, t 1 ) = 0, v(x, T ) = δg δm (x, m(T ))(µ(T )) in T d . ( 10
)
The proof requires several steps. We begin by computing the directional derivative of the map J. For this, let us fix an admissible direction (µ, z). We assume for a while that (µ, z) is smooth and satisfies

∂ t µ -∆µ + div(z) = 0 in T d × (t 0 , T ), µ(t 0 ) = 0. ( 11
)
We also assume that µ = 0 and z = 0 on the time interval [t 0 , t 0 + ε] for some ε > 0. We then consider the map h → J((m, w) + h(µ, z)) for |h| small. As m > 0 on (t 0 , T ] and µ = 0 on [t 0 , t 0 + ε], we have m + hµ > 0 for |h| small enough. Then, by optimality of (m, w) for J, we have:

d dh | h=0 J((m, w) + h(µ, z)) = 0, d 2 dh 2 | h=0 J((m, w) + h(µ, z)) ≥ 0,
where

d dh | h=0 J((m, w) + h(µ, z)) = T t0 T d µL(x, w m ) + mD q L x, w m • ( z m - µw m 2 ) + T t0 T d f (x, m(t))µ + T d g(x, m(T ))µ(T )
and

d 2 dh 2 | h=0 J((m, w) + h(µ, z)) = T t0 T d 2µD q L(x, w m ) • ( z m - µw m 2 ) + mD q L x, w m • (- µz m 2 + 2 µ 2 w m 3 - µz m 2 ) + T t0 T d mD 2 qq L x, w m ( z m - µw m 2 ) • ( z m - µw m 2 ) + T t0 T d T d δf δm (x, m(t), y)µ(x, t)µ(y, t) + T d T d δg δm (x, m(T ), y)µ(x, T )µ(y, T ).
Recalling that w = -mD p H(x, Du), we can rearrange the expressions:

d dh | h=0 J((m, w) + h(µ, z)) = T t0 T d µL(x, w m ) + D q L x, w m • (z + µD p H(x, Du)) + T t0 T d f (x, m(t))µ + T d g(x, m(T ))µ(T ) and d 2 dh 2 | h=0 J((m, w) + h(µ, z)) = T t0 T d m -1 D 2 qq L x, w m (z + µD p H(x, Du)) • (z + µD p H(x, Du)) + T t0 T d T d δf δm (x, m(t), y)µ(x, t)µ(y, t) + T d T d δg δm (x, m(T ), y)µ(x, T )µ(y, T ).
We denote by J (t 0 , m, u; µ, z) the right-hand side of the above expression:

J (t 0 , m, u; µ, z) := T t0 T d m -1 D 2 qq L x, w m (z + µD p H(x, Du)) • (z + µD p H(x, Du)) + T t0 T d T d δf δm (x, m(t), y)µ(x, t)µ(y, t) + T d T d δg δm (x, m(T ), y)µ(x, T )µ(y, T ). Note that J (t 0 , m, u; µ, z) is defined for µ, z ∈ L 2 (T d × (0, T ))
. By regularization, one has therefore: [START_REF] Lasry | Jeux à champ moyen. I. Le cas stationnaire[END_REF] holds in the sense of distribution. Let us also recall [START_REF] Cardaliaguet | The master equation and the convergence problem in mean field games[END_REF] that the map (x, y) → δf δm (x, m, y) is not symmetric, but satisfies the relation:

J (t 0 , m, u; µ, z) ≥ 0 for any µ, z ∈ L 2 (T d × (t 0 , T )) such that
δf δm (x, m, y) = δf δm (y, m, x) + f (x, m) -f (y, m).
So, for any µ, µ ′ such that

T d µ = T d
µ ′ = 0, we have

T d T d δf δm (x, m, y)µ(x)µ ′ (y) = T d T d δf δm (x, m, y)µ(y)µ(x). ( 12
) Lemma 4.4. Let t 1 ∈ [t 0 , T ). For any (µ, z) ∈ L 2 (T d × (t 1 , T )) such that ∂ t µ -∆µ + div(z) = 0 in T d × (t 1 , T ), µ(t 1 ) = 0, ( 13
)
equality J (t 1 , u, m; µ, z) = 0 holds if and only if there exists v ∈ C 0 (T d × (t 1 , T )) such that the pair (v, µ) is a solution to the linearized problem [START_REF] Ladyženskaja | N Linear and quasilinear equations of parabolic type[END_REF] and

z = -µD p H(x, Du) -mD 2 pp H(x, Du)Dv. Proof.
Let us first assume that (v, µ) is a solution to the linearized system [START_REF] Ladyženskaja | N Linear and quasilinear equations of parabolic type[END_REF] and let us set z = -µD p H(x, Du) -mD 2 pp H(x, Du)Dv. Then (µ, z) satisfies ( 13) and

J (t 1 , u, m; µ, z) = T t1 T d mD 2 qq L x, w m D 2 pp H(x, Du)Dv • D 2 pp H(x, Du)Dv + T t1 T d T d δf δm (x, m(t), y)µ(x, t)µ(y, t) + T d T d δg δm (x, m(T ), y)µ(x, T )µ(y, T ).
As H is uniformly convex in the gradient variable and w = -mD p H(x, Du), we have

D 2 qq L x, w m D 2 pp H(x, Du) = I d , (14)
so that

J (t 1 , u, m; µ, z) = T t1 T d mD 2 pp H(x, Du)Dv • Dv + T t1 T d T d δf δm (x, m(t), y)µ(x, t)µ(y, t) + T d T d δg δm (x, m(T ), y)µ(x, T )µ(y, T ) = 0,
where the last equality is obtained by integrating the first equation in [START_REF] Ladyženskaja | N Linear and quasilinear equations of parabolic type[END_REF] multiplied by µ and adding it to the second one multiplied by v.

Conversely let us assume that J (t 1 , u, m; µ, z) = 0 holds for some (µ, z) satisfying ( 13). Let v be the (continuous) solution to

     -∂ t v -∆v + D p H(x, Du) • Dv = δf δm (x, m)(µ) in T d × (t 1 , T ), v(x, T ) = δg δm (x, m(T ))(µ(T )) in T d .
We first claim that (µ, z) is a minimum point of J (t 1 , u, m; •, •) under the constraint [START_REF] Lasry | Mean field games[END_REF]. Indeed, by dynamic programming, (m, w) is optimal for 13) holds, which proves the claim.

J(t 1 , m(t 1 ); •, •). So J (t 1 , u, m; µ ′ , z ′ ) ≥ 0 = J (t 1 , u, m; µ, z) for any (µ ′ , z ′ ) ∈ L 2 such that (
Next we prove that (v, µ) is a solution to the linearized system [START_REF] Ladyženskaja | N Linear and quasilinear equations of parabolic type[END_REF]. As (µ, z) is a minimum point of the quadratic map J (t 1 , u, m; •, •) under the constraint (13), we have by first order necessary condition and for any (µ ′ , z ′ ) such that (13) holds,

T t1 T d m -1 D 2 qq L x, w m (z + µD p H(x, Du)) • (z ′ + µ ′ D p H(x, Du)) + T t1 T d T d δf δm (x, m(t), y)µ(x, t)µ ′ (y, t) + T d T d δg δm (x, m(T ), y)µ(x, T )µ ′ (y, T ) = 0.
(we used [START_REF] Lasry | Jeux à champ moyen. II. Horizon fini et contrôle optimal[END_REF] in the above equality). Computing as usual d dt T d vµ ′ and integrating in time we find:

T d µ ′ (T ) δg δm (x, m(T ))(µ(T )) = T t1 T d µ ′ (D p H(x, Du) • Dv - δf δm (x, m)(µ)) + z ′ • Dv.
Recalling [START_REF] Lasry | Jeux à champ moyen. II. Horizon fini et contrôle optimal[END_REF], this implies that

T t1 T d m -1 D 2 qq L x, w m (z + µD p H(x, Du)) • (z ′ + µ ′ D p H(x, Du)) + T t1 T d (z ′ + µ ′ D p H(x, Du)) • Dv = 0.
Let us note that, for any α ∈ L 2 , there exists a unique solution µ ′ ∈ L 2 to

∂ t µ ′ -∆µ ′ -div(µ ′ D p H(x, Du)) + div(α) = 0 in T d × (t 1 , T ), µ(t 1 ) = 0.
14 Then the pair (µ ′ , z ′ ), with z ′ := -µ ′ D p H(x, Du) + α, satisfies [START_REF] Lasry | Mean field games[END_REF], so that we have

T t1 T d m -1 D 2 qq L x, w m (z + µD p H(x, Du)) + Dv • α = 0 for any α ∈ L 2 . Therefore m -1 D 2 qq L x, w m (z + µD p H(x, Du)) + Dv = 0,
that, thanks to [START_REF] Lions | Sur l'unicité rétrograde dans les problèmes mixtes paraboliques[END_REF], can be rewritten as

z = -µD p H(x, Du) -mD 2 pp H(x, Du)Dv.
This shows that the pair (v, µ) is a solution to the linearized system [START_REF] Ladyženskaja | N Linear and quasilinear equations of parabolic type[END_REF].

Proof of Theorem 4.3. Let (v, µ) be a solution of the linearized system [START_REF] Ladyženskaja | N Linear and quasilinear equations of parabolic type[END_REF]. From Lemma 4.4, one has [START_REF] Lasry | Jeux à champ moyen. I. Le cas stationnaire[END_REF] and J (t 0 , u, m; µ, z) = 0. Using again Lemma 4.4, we know that there exists ṽ such that (ṽ, µ) is a solution to the linearized system (10) on T d × (t 0 , T ) with initial condition µ(t 0 ) = 0 and z = -µD p H(x, Du)-mD 2 pp H(x, Du)Dṽ. By parabolic regularity, as (ṽ, µ) solves [START_REF] Ladyženskaja | N Linear and quasilinear equations of parabolic type[END_REF], it satisfies ṽ ∈ C 1 and µ ∈ C 0 . Thus z(t 1 ) = 0 = -µ(t 1 )D p H(x, Du(t 1 )) -m(t 1 )D 2 pp H(x, Du(t 1 ))Dṽ(t 1 ). Now recalling that µ(t 1 ) = 0, m(t 1 ) > 0 and D 2 pp H > 0, this implies that Dṽ(t 1 ) = 0. So the pair (Dṽ, µ) is a classical solution to

J (t 1 , u, m; µ, z) = 0, where z = -µD p H(x, Du) -mD 2 pp H(x, Du)Dv. Let us extend (µ, z) to T d × [t 0 , T ] by setting (µ, z) = 0 on T d × [t 0 , t 1 ]. Then (µ, z) satisfies
   -∂ t (∂ xi ṽ) -∆(∂ xi ṽ) + g i = 0 in T d × (t 1 , T ), i = 1, . . . , d, ∂ t µ -∆µ + div(h) = 0 in T d × (t 1 , T ) µ(x, t 1 ) = 0, Dṽ(x, t 1 ) = 0 in T d .
where

g i (x, t) = ∂ xi D p H(x, Du) • Dṽ - δf δm (x, m)(µ) and h(x, t) = -µD p H(x, Du) -mD 2 pp H(x, Du)Dṽ. Note that d i=1 |g i (x, t)| 2 ≤ C |Dṽ| 2 + |D 2 ṽ| 2 + µ(t) 2 L 2 while |h(x, t)| 2 ≤ C |µ(x, t)| 2 + |Dṽ(x, t)| 2 and |div(h)(x, t)| 2 ≤ C |Dṽ| 2 + |D 2 ṽ| 2 + |µ(x, t)| 2 + |Dµ(x, t)| 2 .
Then Theorem 6.1 in the Appendix states that Dṽ = 0 and µ = 0. As (v, µ) solves [START_REF] Ladyženskaja | N Linear and quasilinear equations of parabolic type[END_REF], this also implies that v = 0.

Application to a learning procedure

In order to illustrate the notion of stable solution in potential games, we consider a learning procedure and show that this procedure converges if one starts from a neighborhood of a stable MFG equilibrium. In [START_REF] Cardaliaguet | Learning in Mean Field Games: the Fictitious Play[END_REF] the following algorithm (inspired by the Fictitious Play [START_REF] Brown | Iterative solution of games by Fictitious Play[END_REF]) was introduced: let µ 0 ∈ C 0 ([0, T ], P(T d )) and define by induction (u n , m n ) the solution to

   -∂ t u n+1 -∆u n+1 + H(x, Du n+1 ) = f (x, µ n )) ∂ t m n+1 -∆m n+1 -div(m n+1 D P H(x, Du n+1 )) = 0 m n+1 (0) = m 0 , u n+1 (T, x) = g(x, µ n (T )) (15)
and

µ n+1 = n n + 1 µ n + 1 n + 1 m n+1 .
Following ideas of Monderer and Shapley [START_REF] Monderer | Fictitious play property for games with identical interests[END_REF], it was proved in [START_REF] Cardaliaguet | Learning in Mean Field Games: the Fictitious Play[END_REF] that the family (u n , m n ) is compact in C 1,0 × C 0 and that any converging subsequence is a solution of the MFG system.

Our main result is the convergence of (u n , m n ) to a stable solution (u, m) as soon as µ 0 is sufficiently close to m. This can be understood as a robustness property of the equilibrium (u, m).

Theorem 5.1. Assume that the game is potential, i.e., satisfies (5) for some C 1 maps F, G : P(T d ) → R. Let (u, m) be a solution to the MFG system (1) starting from (0, m 0 ) with m 0 ∈ P(T d ) and assume that the solution is stable. Then there exists δ > 0 such that, if sup

t∈[0,T ] d 1 (m(t), µ 0 (t)) ≤ δ, the sequence (u n , m n ) defined above converges to (u, m) as n → +∞ in C 1,0 × C 0 .
In other worlds, stable equilibria are local attractors for the learning procedure. Note also that, in contrast to [START_REF] Cardaliaguet | Learning in Mean Field Games: the Fictitious Play[END_REF], we show here that the full sequence (u n , m n ) converges, although the MFG system may have several solutions.

We now start the proof of Theorem 5.1. Let us recall that, by the stability assumption of the MFG equilibrium (u, m), the linearized system

   -∂ t v -∆v + D p H(x, Du) • Dv = δf δm (x, m)(ρ) in T d × (0, T ), ∂ t ρ -∆ρ -div(ρD p H(x, Du)) -div(mD 2 pp H(x, Du)Dv) = 0 in T d × (0, T ) ρ(x, 0) = 0, v(x, T ) = δg δm (x, m(T ))(ρ(T )) in T d . (16)
has the pair (v, ρ) = (0, 0) as unique solution. As a consequence, we have:

Lemma 5.2. There is a constant C > 0 such that, for any a, b ∈ C 0 (T d × [0, T ]), c ∈ C 0 (T d ) and (v, ρ) solution to    -∂ t v -∆v + D p H(x, Du) • Dv = δf δm (x, m)(ρ) + a(x, t) in T d × (0, T ), ∂ t ρ -∆ρ -div(ρD p H(x, Du)) -div(mD 2 pp H(x, Du)Dv) = div(b(x, t)) in T d × (0, T ) ρ(x, 0) = 0, v(x, T ) = δg δm (x, m(T ))(ρ(T )) + c in T d , (17)
one has v C 1,0 + ρ C 0 ≤ C ( a C 0 + b C 0 + c C 0 ) .
Proof. We argue by contradiction. Since the system is affine, this means that there exists a n , b n , c n and (v n , ρ n ) solution to (17) with

a n C 0 + b n C 0 + c n C 0 ≤ 1 and θ n := v n C 1,0 + ρ n C 0 ≥ n.
Let us set ṽn := v n θ n , ρn := ρ n θ n . The pair (ṽ n , ρn ) solves

   -∂ t ṽn -∆ṽ n + D p H(x, Du) • Dṽ n = δf δm (x, m)(ρ n ) + a n (x, t)/θ n in T d × (0, T ), ∂ t ρn -∆ρ n -div(ρ n D p H(x, Du)) -div(mD 2 pp H(x, Du)Dṽ n ) = div (b n (x, t)/θ n ) in T d × (0, T ) ρn (x, 0) = 0, ṽn (x, T ) = δg δm (x, m(T ))(ρ(T )) + c n /θ n in T d .
The ṽn satisfy a linear parabolic equation with bounded coefficients. Therefore (see for instance, Theorem 9.1 and Theorem 10.1 in [START_REF] Ladyženskaja | N Linear and quasilinear equations of parabolic type[END_REF]), the sequences (ṽ n ) and (Dṽ n ) are bounded in C α,α/2 for some α ∈ (0, 1). In the same way, the ρn satisfy a linear parabolic equation in divergence form with bounded coefficients, so that the sequence (ρ n ) is bounded in C α,α/2 . Thus, up to a subsequence (ṽ n ) and (Dṽ n ) converge uniformly to some v and Dv while (ρ n ) converges to some ρ. Note that, by definition of θ n and (ṽ n , ρn ), we have

v C 1,0 + ρ C 0 = lim n→+∞ ( ṽn C 1,0 + ρn C 0 ) = 1.
On the other hand, the pair (v, ρ) is a weak solution to the linearized system (16), so that, by our assumption on this system, (v, ρ) must actually vanish. So there is a contradiction.

Lemma 5.3. There exists η > 0 and C > 0 such that, if

D(u n+1 -u) C 0 + m n+1 -m C 0 ≤ η, then u n+1 -u C 1,0 + m n+1 -m C 0 ≤ C m n+1 -µ n C 0 .
Proof. We rewrite system [START_REF] Monderer | Fictitious play property for games with identical interests[END_REF] as (17

) with v = u n+1 -u, ρ = m n+1 -m, a := H(x, Du n+1 ) -H(x, Du) -D p H(x, Du) • D(u n+1 -u) +f (x, µ n (t)) -f (x, m(t)) - δf δm (m(t))(m n+1 (t) -m(t)), b := m n+1 D p H(x, Du n+1 ) -mD p H(x, Du) -(m n+1 -m)D p H(x, Du) -mD 2 pp H(x, Du)D(u n+1 -u) and c(x) := g(x, µ n (T )) -g(x, m(T )) - δg δm (m(t))(m n+1 (T ) -m(T )).
Then, by Lemma 5.2, we get

u n+1 -u C 1,0 + m n+1 -m C 0 ≤ C( a C 0 + b C 0 + c C 0 ) ≤ C D(u n+1 -u) 2 C 0 + µ n -m 2 C 0 + m n+1 -µ n C 0 + m n+1 -m 2 C 0 ≤ C D(u n+1 -u) 2 C 0 + m n+1 -m 2 C 0 + m n+1 -µ n C 0
where the constant C in the last inequality depends only on the data. Thus, if η > 0 is small enough and D(u n+1 -u) C 0 + m n+1 -m C 0 is not larger than η, one can absorb the square terms of the right-hand side into the left-hand side to get the conclusion.

Proof of Theorem 5.1. Following [START_REF] Cardaliaguet | Learning in Mean Field Games: the Fictitious Play[END_REF], the difference (m n+1 -µ n ) converges to 0 in C 0 : so we can find N ≥ 0 large enough so that

m n+1 -µ n C 0 ≤ η/(2C) ∀n ≥ N, (18)
where η and C are given by Lemma 5.3.

Next we note that the u n are uniformly Lipschitz continuous. So the m n solve a parabolic equation in divergence form with bounded coefficient and therefore are bounded in C α,α/2 . Thus so are the µ n . Moreover, by definition of µ n , we have that

µ n+1 -µ n C 0 = 1 n m n+1 -µ n C 0 ≤ C n .
As u n+2 and u n+1 solve the same equation with right-hand side given respectively by f (x, µ n ) and f (x, µ n+1 ) and terminal condition given by g(x, µ n (T )) and g(x, µ n+1 (T )) respectively, we conclude that

u n+2 -u n+1 C 1,0 ≤ C f (•, µ n ) -f (•, µ n+1 ) C 1,0 + g(•, µ n (T )) -g(•, µ n+1 (T )) C 1,0 ≤ C n .
(19) Plugging this estimate into the equation satisfied by m n+1 -m n+2 also gives

m n+1 -m n+2 C 0 ≤ C n . (20)
From now on we assume that N is so large that the right-hand sides C/n in ( 19) and (20) are less than η/4 for any n ≥ N . Next we note that, if µ 0 is sufficiently close to m, then u 1 is close to u (in C 1,0 ) and m 1 is close to m (in C 0 ). This in turn shows that µ 1 is close to m (in C 0 ). By induction, we obtain that, for any fixed n, one can choose µ 0 sufficiently close to m so that u n+1 is close to u (in C 1,0 ) and m n+1 is close to m (in C 0 ). Applying this argument to the integer N defined above, we obtain that one can choose µ 0 so close to m that

u N +1 -u C 1,0 + m N +1 -m C 0 ≤ η.
We claim that this inequality propagate to any n + 1 ≥ N + 1. Indeed, let us assume that it holds for some n + 1. Then, by Lemma 5.3, we have

u n+1 -u C 1,0 + m n+1 -m C 0 ≤ C m n+1 -µ n C 0 ≤ η/2
thanks to (18). So, by ( 19) and (20) and the choice of N , we also have

u n+2 -u C 1,0 + m n+2 -m C 0 ≤ u n+1 -u C 1,0 + m n+1 -m C 0 + u n+2 -u n+1 C 1,0 + m n+1 -m n+2 C 0 ≤ η,
which proves the claim.

In particular, we can apply Lemma 5.3 to any n ≥ N to get

u n+1 -u C 1,0 + m n+1 -m C 0 ≤ C m n+1 -µ n C 0 ,
where the right-hand side tends to 0 as n → +∞ by Theorem 2.1 in [START_REF] Cardaliaguet | Learning in Mean Field Games: the Fictitious Play[END_REF]. This proves the convergence of (u n ) and (m n ) to u and m respectively.

Appendix

Let (u, µ) = ((u i ) i=1,...,d , µ) be a classical solutions to the forward-backward system

-∂ t u i -∆u i + g i (x, t) = 0 in T d × (0, T ), i = 1, . . . , d, ∂ t µ -∆µ + div(h(x, t)) = 0 in T d × (0, T ) We suppose that u i (x, 0) = µ(x, 0) = 0 in T d , i = 1, . . . , d, d i=1 |g i (x, t)| 2 ≤ C 0 µ(•, t) 2 2 + |u(x, t)| 2 + |Du(x, t)| 2 , |h(x, t)| 2 ≤ C 0 |µ(x, t)| 2 + |u(x, t)| 2 , |div(h)(x, t)| 2 ≤ C 0 |µ(x, t)| 2 + |Dµ(x, t)| 2 + |u(x, t)| 2 + |Du(x, t)| 2 .
Theorem 6.1. Under the above assumptions, one has u(x, t) = 0 and µ(x, t) = 0 on

T d × [0, T ].
The argument is standard and goes back to J.L. Lions and B. Malgrange in [START_REF] Lions | Sur l'unicité rétrograde dans les problèmes mixtes paraboliques[END_REF] (see also Cannarsa-Tessitore [START_REF] Cannarsa | Optimality conditions for boundary control problems of parabolic type[END_REF] for a forward-backward system).

Proof. It is enough to argue for T > 0 sufficiently small and prove that µ = 0 and u = 0 on [0, T /2]. Let θ : [0, T ] → [0, 1] be a smooth, non increasing function with θ

(t) = 1 in [0, T /2], θ(t) = 0 in [2T /3, T ] and θ ′ (t) ∞ ≤ C/T . For κ ≥ 1 we set ũi (x, t) = e κ(t-T ) 2 /2 θ(t)u i (x, t), μ(x, t) = e κ(t-T ) 2 /2 θ(t)µ(x, t).
Then (ũ, μ) = ((ũ i ), μ) satisfies    (i) -∂ t ũi -∆ũ i + κ(t -T )ũ i + e κ(t-T ) 2 /2 θ ′ u i + e κ(t-T ) 2 /2 θg i = 0, (ii) ∂ t μ -∆μ -κ(t -T )μ -e κ(t-T ) 2 /2 θ ′ µ + e κ(t-T ) 2 /2 θdiv(h) = 0, (iii) ũi (x, 0) = μ(x, 0) = ũi (x, T ) = μ(x, T ) = 0. 21)-(ii) by ∂ t μ and integrating in time-space, we obtain, after integration by parts the second term:

T 0 T d (∂ t μ) 2 + 1 2 ∂ t |D μ| 2 - κ 2
(t -T )∂ t (μ) 2 -∂ t μ e κ(t-T ) 2 /2 θ ′ µ -e κ(t-T ) 2 /2 θdiv(h) = 0.

We integrate the second term in time and integrate by parts in time the third one: taking into account the fact that μ(•, 0) = μ(•, T ) = 0, we get:

T 0 T d (∂ t μ) 2 + κ 2 (
μ) 2 -∂ t μ e κ(t-T ) 2 /2 θ ′ µ -e κ(t-T ) 2 /2 θdiv(h) = 0. By Young's inequality we obtain

T 0 T d 1 2 (∂ t μ) 2 + κ 2 (μ) 2 ≤ C T 0 T d e κ(t-T ) 2 (θ ′ ) 2 µ 2 + e κ(t-T ) 2 θ 2 |div(h)| 2 . ( 22
)
We argue in the same way for u i : we multiply (21)-(i) by ∂ t ũi and integrate in space-time:

T 0 T d -(∂ t ũi ) 2 + 1 2 ∂ t |D ũi | 2 + κ 2
(t -T )∂ t (ũ i ) 2 + ∂ t ũi e κ(t-T ) 2 /2 θ ′ u i + e κ(t-T ) 2 /2 θg i = 0.

This yields to

T 0 T d 1 2 (∂ t ũi ) 2 + κ 2 (ũ i ) 2 ≤ C T 0 T d
e κ(t-T ) 2 (θ ′ ) 2 u 2 i + e κ(t-T ) 2 θ 2 g 2 i .

(23)

Next we multiply (21)-(ii) by μ and integrate in time-space: using as usual that μ(•, 0) = μ(•, T ) = 0, we obtain

T 0 T d
|D μ| 2 -κ(t -T )(μ) 2 + μ(-e κ(t-T ) 2 /2 θ ′ µ) -D μ • (e κ(t-T ) 2 /2 θh) = 0.

Hence, by Young's inequality,

T 0 T d 1 2 |D μ| 2 ≤ T 0 T d μ2 + C e κ(t-T ) 2 (θ ′ ) 2 µ 2 + e κ(t-T ) 2 θ 2 h 2 . ( 24
)
In the same way, we multiply ( 21)-(i) by u i and integrate in time-space:

T 0 T d |D ũi | 2 + κ(t -T )(ũ i ) 2
+ ũi e κ(t-T ) 2 /2 θ ′ u i + e κ(t-T ) 2 /2 θg i = 0.

Then, for ε ∈ (0, 1) to be chosen later,

T 0 T d |D ũi | 2 ≤ T 0 T d (κT + ε -1 )(ũ i ) 2 + Cε e κ(t-T ) 2 (θ ′ ) 2 u 2 i + e κ(t-T ) 2 θ 2 g 2 i . (25)
We now take into account the assumptions on the g i , h and div(h): (22) becomes

κ 2 T 0 μ(t) 2 2 ≤ C
T 0 e κ(t-T ) 2 (θ ′ ) 2 µ(t) 2 2 + e κ(t-T ) 2 θ 2 u(t) We can now fix T > 0 small enough, so that, for any κ large enough, Dividing by e κ(T /2) 2 and letting κ → +∞ yields to µ = 0 and u = 0 on [0, T /2].

κ 2 T 0 μ(t)

  ) over i yields, for ε small enough (depending only on the constant C in (29)), + e κ(t-T ) 2 (θ ′ ) 2 u(t) 2 2 + μ(t) 2 H 1 .

	Summing (29T 0 2 Plugging (28) into the above inequality gives: D ũ(t) 2 2 ≤ C T 0 (κT + 1) ũ(t) 2
			0	T	D ũ(t) 2 2 ≤ C	0	T	(κT + 1) ũ(t) 2 2 + μ(t) 2 2 + e κ(t-T ) 2 (θ ′ ) 2 ( u(t) 2 2 + µ(t) 2 2 ).	(30)
		Collecting (26), (27) (summing on i = 1, . . . , d), (28) and (30) yields to
	κ	0	T	μ(t) 2 2 + ũ(t) 2 2 ≤ C ≤ C	0 0	T T	e κ(t-T ) 2 (θ ′ ) 2 ( µ(t) 2 2 + u(t) 2 2 ) + ũ(t) 2 H 1 + μ(t) 2 H 1 e κ(t-T ) 2 (θ ′ ) 2 ( µ(t) 2 2 + u(t) 2 2 ) + (κT + 1) ũ(t) 2 2 + μ(t) 2 2
										2 2 + Du(t) 2 2 + µ(t) 2 2 + Dµ(t) 2 2	.
	Rearranging we obtain
					κ	0	T	μ(t) 2 2 ≤ C	0	T	e κ(t-T ) 2 (θ ′ ) 2 µ(t) 2 2 + ũ(t) 2 H 1 + μ(t) 2 H 1 .	(26)
	Note that we used the notation v(t) 2 H 1 := v(t) 2 2 + Dv(t) 2 2 . The same argument for (23)
	yields to	κ	0	T	ũi (t) 2 2 ≤ C	0	T	e κ(t-T ) 2 (θ ′ ) 2 u i (t) 2 2 + ũ(t) 2 H 1 + μ(t) 2 H 1 ,	(27)
	while (24) and (25) become respectively:
							0	T	D μ(t) 2 2 ≤ C	0	T	e κ(t-T ) 2 (θ ′ ) 2 µ(t) 2 2 + ũ(t) 2 2 + μ(t) 2 2 ,	(28)
	and
	0	T		D ũi (t) 2 2 ≤		0	T	(κT + ε -1 ) ũi (t) 2 2 + Cε e κ(t-T ) 2 (θ ′ ) 2 u i (t) 2 2 + ũ(t) 2 H 1 + μ(t) 2 H 1 . (29)
										20