Jules Berry
email: jules.berry@insa-rennes.fr

Olivier Ley
email: olivier.ley@insa-rennes.fr

Francisco J Silva
email: francisco.silva@unilim.fr

Approximation and perturbations of stable solutions to a stationary mean field game system

Keywords: 2020 Mathematics subject classification. 35J47, 35Q89, 65N20 Mean field games, Newton method, finite element method, numerical methods, non-uniqueness of solutions

This work introduces a new general approach for the numerical analysis of stable equilibria to second order mean field games systems in cases where the uniqueness of solutions may fail. For the sake of simplicity, we focus on a simple stationary case. We propose an abstract framework to study these solutions by reformulating the mean field game system as an abstract equation in a Banach space. In this context, stable equilibria turn out to be regular solutions to this equation, meaning that the linearized system is well-posed. We provide three applications of this property: we study the sensitivity analysis of stable solutions, establish error estimates for their finite element approximations, and prove the local converge of Newton's method in infinite dimensions.

Introduction

Mean Field Games (MFG for short) were introduced independently by Lasry-Lions [START_REF] Lasry | Mean Field Games[END_REF][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] and Huang-Caines-Malhamé [START_REF] Huang | Large Population Stochastic Dynamic Games: Closed-Loop McKean-Vlasov Systems and the Nash Certainty Equivalence Principle[END_REF][START_REF] Huang | Large-Population Cost-Coupled LQG Problems With Nonuniform Agents: Individual-Mass Behavior and Decentralized ε-Nash Equilibria[END_REF]. The goal of this theory is to study (stochastic) differential games with a large number of interchangeable players. We refer the reader to [START_REF] Achdou | Mean Field Games: Cetraro[END_REF][START_REF] Carmona | Probabilistic Theory of Mean Field Games with Applications I-II. Probability Theory and Stochastic Modelling[END_REF][START_REF] Bensoussan | Mean Field Games and Mean Field Type Control Theory[END_REF][START_REF] Gomes | Mean Field Games Models-A Brief Survey[END_REF][START_REF] Lions | Cours Au College de France[END_REF] for general references on this topic.

The numerical analysis of MFG systems introduced in [START_REF] Lasry | Mean Field Games[END_REF][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] has been extensively studied under a monotonicity assumption also introduced by Lasry and Lions, see [START_REF] Achdou | Mean Field Games: Cetraro[END_REF]Chapter 4], [START_REF] Achdou | Finite Difference Methods for Mean Field Games[END_REF], and the references therein. Indeed, the latter provides a sufficient condition for the uniqueness of solutions to MFG systems, which allows to show the convergence of numerical methods in [START_REF] Achdou | Mean Field Games: Convergence of a Finite Difference Method[END_REF] and error estimates in [START_REF] Bonnans | Error Estimates of a Theta-Scheme for Second-Order Mean Field Games[END_REF]. In the absence of this monotonicity assumption, uniqueness may fail (see [START_REF] Briani | Stable Solutions in Potential Mean Field Game Systems[END_REF][START_REF] Bardi | On non-uniqueness and uniqueness of solutions in finite-horizon mean field games[END_REF][START_REF] Cirant | Time-dependent focusing mean-field games: the sub-critical case[END_REF][START_REF] Cirant | On the existence of oscillating solutions in non-monotone mean-field games[END_REF]) and the study of MFG systems with several solutions is delicate both from the theoretical and numerical points of view. In [START_REF] Briani | Stable Solutions in Potential Mean Field Game Systems[END_REF], Briani and Cardaliaguet define a particular notion of solution for second order potential MFG systems, the so-called stable solutions (see also [START_REF] Bressan | Generic properties of first-order mean field games[END_REF] for a related notion in the context of first order mean field games). These solutions may not be unique but the authors show in [START_REF] Briani | Stable Solutions in Potential Mean Field Game Systems[END_REF] that they have some interesting properties motivating their name: stable solutions are isolated and the fictitious play algorithm, introduced in Cardaliaguet-Hadikhanloo [START_REF] Cardaliaguet | Learning in Mean Field Games: The Fictitious Play[END_REF], converges locally to these solutions. We also mention the recent work by Tang-Song [START_REF] Tang | Learning Optimal Policies in Potential Mean Field Games: Smoothed Policy Iteration Algorithms[END_REF], where the authors implement a smoothed policy iteration method to locally approximate stable solutions. In this paper, we provide new results in this direction, which reinforce the importance of the notion of stable solutions. We are going to prove that stable solutions are indeed stable under perturbations and that local convergence holds for their approximations by finite element methods and Newton iterations. This paper is the first in a series of works dealing with the numerical analysis of stable equilibria to MFG models. Our goal is to introduce a general framework that covers different types of MFG systems under fairly general assumptions. In order to convey our main ideas, we focus in this paper on the following stationary MFG system

-∆u + 1 2 |Du| 2 + λu = f (m) in T d , -∆m -div (mDu) + λm = λm 0 in T d , (1)
where λ > 0 is a given constant and m 0 : T d → R and f : R → R are given functions. This system has been introduced in the monograph by Bensoussan-Frehse-Yam [5, Chapter 7] and has been furtherly studied in [START_REF] Ferreira | Existence of weak solutions to stationary mean-field games through variational inequalities[END_REF][START_REF] Gomes | The selection problem for some first-order stationary mean-field games[END_REF]. We briefly recall its interpretation in Section 2.

In this work we give a definition of stable solutions for system (1) (Definition 2.5 below) and prove that, except for at most finitely many values of λ > 0, every solution to (1) is stable. Furthermore, if the coupling f is monotone, then the unique solution to (1) is stable for every λ > 0. As in [START_REF] Briani | Stable Solutions in Potential Mean Field Game Systems[END_REF], we prove that stable solutions are isolated and we are able to slightly improve the result in [START_REF] Briani | Stable Solutions in Potential Mean Field Game Systems[END_REF] by considering weaker norms.

Let us now present the main contribution of this paper. We reformulate system [START_REF] Achdou | Finite Difference Methods for Mean Field Games[END_REF] in the form

F (u, m) = 0,
where F : X → X is a nonlinear mapping, defined on a Banach space X, having the form

F = I + T • G.
More precisely, given a suitable Banach space Z, we choose T : Z → X as the linear operator defined by T (f, g) = (v, ρ), where (v, ρ) solves

-∆v + λv = f in T d , -∆ρ + λρ = g in T d ,
and G : X → Z the nonlinear mapping defined by

G(v, ρ) = 1 2 |Dv| 2 -f (ρ), -div (ρDu) -λm 0 .
We refer the reader to Section 3 for the details of this reformulation. In the case where the mapping G is differentiable, the mapping F is also differentiable and the stability of a solution (u, m) to ( 1) is equivalent to the injectivity of the differential

dF [u, m] = I + T • dG[u, m].
Thus, if T is a compact operator and (u, m) is a stable solution to [START_REF] Achdou | Finite Difference Methods for Mean Field Games[END_REF], then dF [u, m] is an injective perturbation of the identity by a compact linear operator. Henceforth, by the Fredholm alternative, we deduce that dF [u, m] is an isomorphism on X. This isomorphism property will be rigorously established below for Banach spaces of the form

X = C 2,γ (T d ) × C 2,β (T d ) and X = W 1,p (T d ) × L q (T d ),
but many other choices are possible depending on the application in mind.

We now describe three applications of the above isomorphism property for stable solutions.

A first and straightforward application, which follows from the implicit function theorem, concerns the sensitivity analysis of stable solutions to [START_REF] Achdou | Finite Difference Methods for Mean Field Games[END_REF] under perturbations of the coupling function f and the distribution m 0 (Proposition 5.1 below).

In the second application, we make use of the Brezzi-Rappaz-Raviart theory on the approximation of nonlinear problems (see [START_REF] Brezzi | Finite Dimensional Approximation of Nonlinear Problems: Part I: Branches of Nonsingular Solutions[END_REF][START_REF] Girault | Finite Element Methods for Navier-Stokes Equations[END_REF][START_REF] Crouzeix | On numerical approximation in bifurcation theory[END_REF][START_REF] Caloz | Numerical analysis for nonlinear and bifurcation problems[END_REF]) to obtain existence and error estimates for finite element approximations of stable solutions (Theorem 5.4 below). The finite element approximation of a MFG system similar to [START_REF] Achdou | Finite Difference Methods for Mean Field Games[END_REF] has been studied by Osborne and Smears in [START_REF] Osborne | Analysis and Numerical Approximation of Stationary Second-Order Mean Field Game Partial Differential Inclusions[END_REF] (see [START_REF] Osborne | Finite element approximation of time-dependent Mean Field Games with non-differentiable Hamiltonians[END_REF] for a parabolic MFG system), where the convergence of the approximations is established by using compactness arguments which do not provide error estimates. While finishing this paper, we have learnt about the recent work [START_REF] Osborne | Near and full quasi-optimality of finite element approximations of stationary second-order mean field games[END_REF] by the same authors addressing this issue. The results in [START_REF] Osborne | Analysis and Numerical Approximation of Stationary Second-Order Mean Field Game Partial Differential Inclusions[END_REF][START_REF] Osborne | Near and full quasi-optimality of finite element approximations of stationary second-order mean field games[END_REF] deal only with the case where the coupling term f is (strongly) monotone, which ensures the uniqueness of the solution to the MFG system. In contrast, our results apply locally around any stable solution without requiring any monotonicity of the coupling and rely on a completely different approach.

In our last application, we provide convergence rates for the iterates of Newton's method in infinite dimension applied to (1) in various functional spaces (Theorems 5.8, 5.10 and 5.11 below). Let us also mention that the analysis of Newton's method in infinite dimensions to approximate the solution to a time-dependent MFG system, with monotone couplings, has been recently addressed by Camilli and Tang [START_REF] Camilli | A Convergence Rate for the Newton's Method for Mean Field Games with Non-Separable Hamiltonians[END_REF] using different techniques. Compared with their approach, our result follows directly from classical convergence results of Newton's iterates in function spaces (see [START_REF] Hinze | Optimization with PDE Constraints[END_REF][START_REF] Zeidler | Nonlinear Functional Analysis and Its Applications: Fixed-point Theorems[END_REF][START_REF] Dontchev | Implicit Functions and Solution Mappings: A View from Variational Analysis[END_REF]) which allow us to deal with non-monotone couplings and to obtain convergence rates in stronger norms.

The paper is structured as follows. In Section 2 we study the well-posedness of (1) and provide some useful estimates on the solutions. We then introduce the definition of stable solution to [START_REF] Achdou | Finite Difference Methods for Mean Field Games[END_REF], which is similar to the one proposed in [START_REF] Briani | Stable Solutions in Potential Mean Field Game Systems[END_REF], and give sufficient conditions for the existence of such solutions. In Section 3 we reformulate solutions to (1) as zeros of a well-chosen nonlinear mapping F and prove the isomorphism property of the differential of F at stable solutions under the assumption that F is continuously differentiable. This assumption is rigorously justified in Section 4 in the case of Hölder and Sobolev spaces. Finally, applications to the sensitivity analysis of (1), to its finite element approximation, and to the convergence of Newton's method, are studied in Section 5.

Notations For k ∈ N and α ∈ (0, 1], we write C k,α (T d ) the usual Hölder space on T d , i.e.

C k,α (T d ) = u ∈ C k (T d ) : u C k,α < +∞ , where u C k,α = |j|≤k ∂ j u L ∞ + |j|=k ∂ j u α , with [u] α = sup y,x∈T d x =y |u(x) -u(y)| |x -y| α .
Similarly, C k,α loc (R) is the space of locally Hölder continuous functions on R, i.e.

f ∈ C k,α loc (R) if f ∈ C k,α (Ω)
for every bounded open subset Ω ⊂ R. We also write C k b (R) the set of k times continuously differentiable functions f on R such that the derivatives f (j) , for 0 ≤ j ≤ k, are bounded. We denote by P(T d ) the space on probability measures over T d and we always identify a measure m ∈ P(T d ) with its density, which we also denote m, provided that the latter exists. For 1 < p ≤ ∞ the dual of the Sobolev space W 1,p (T d ) is denoted by W -1,p ′ (T d ), where 1/p + 1/p ′ = 1, and we reserve the notation H -1 (T d ) for the dual of H 1 (T d ).

For Banach spaces X and Y and a mapping Φ : X → Y , we write dΦ[x] the Fréchet differential of Φ at x ∈ X, when it exists. We also use the notation Y ֒→ X when Y is continuously embedded in X.

By convention, we do not specify integration domains when integrals are considered on T d , i.e.,

f dx := T d f dx.
2 The mean field game system

In this section we establish some properties of the MFG system (1). We first state existence of solutions to (1) as well as a uniqueness result under the Lasry-Lions monotonicity condition on the coupling f . We then define stable solutions to (1) following [START_REF] Briani | Stable Solutions in Potential Mean Field Game Systems[END_REF] and prove that they are isolated. Finally, we provide sufficient conditions to ensure that any classical solution to (1) is stable.

Let us begin by describing the mean field game interpretation of system [START_REF] Achdou | Finite Difference Methods for Mean Field Games[END_REF]. We consider a typical player whose dynamics is governed by the following controlled stochastic differential equation

dX α t = α(X α t )dt + √ 2dB t for all t > 0, X α 0 = x,
where B t is a d-dimensional Brownian motion and α is a feedback control. Assume that the player forecasts ρ as being the evolution of the distribution of players, from which it is possible to compute the weighted averaged density m defined by

m(x) = λ +∞ 0 ρ(t, x)e -λt dt (2)
for some given λ > 0. Then this player aims to minimize the following cost

J(x, α) = E +∞ 0 |α(X α t )| 2 2 + f ( m(X α t )) e -λt dt .
This yields the Hamilton-Jacobi equation for the value function u(x) := inf α J(x, α) :

-∆u + |Du| 2 2 + λu = f ( m).
Since all the players are assumed to be interchangeable, they should all play according to the optimal strategy provided by the Hamilton-Jacobi equation, i.e. α ⋆ (x) = -Du(x) (see [START_REF] Fleming | Deterministic and Stochastic Optimal Control[END_REF][START_REF] Fleming | Controlled Markov Processes and Viscosity Solutions[END_REF][START_REF] Yong | Stochastic Controls[END_REF]). The player must then update the forecasted density by solving the Fokker-Planck equation

∂ t ρ -∆ρ -div(ρDu) = 0, ρ(0) = m 0 , (3)
where m 0 is the initial distribution of agents. This yields and updated averaged density m through [START_REF] Achdou | Mean Field Games: Convergence of a Finite Difference Method[END_REF]. Using integration by parts in (2) one can easily derive the following equation on m -∆m -div(mDu) + λm = λm 0 .

As usual we define an MFG equilibrium as being a fix point of this procedure which corresponds to a solution to [START_REF] Achdou | Finite Difference Methods for Mean Field Games[END_REF]. Notice that once we have a solutions to (1) it is possible to recover the probability density function ρ by solving the Fokker-Planck equation [START_REF] Achdou | Mean Field Games: Cetraro[END_REF].

In all of this paper we make the assumption that m 0 ∈ P(T d ) ∩ C 0,α (T d ) for some α ∈ (0, 1).

(5)

Well-posedness

The following result was established in [START_REF] Bensoussan | Mean Field Games and Mean Field Type Control Theory[END_REF]Chapter 7] by approximation. For the sake of self-containedness we provide a proof in Appendix A based on the classical method introduced by Lasry-Lions.

Theorem 2.1. Assume that f ∈ W 1,∞ (R). Then there exists a classical solution

(u, m) ∈ C 2,α (T d ) × C 2,α (T d ) to (1)
, where the constant α is fixed in [START_REF] Bensoussan | Mean Field Games and Mean Field Type Control Theory[END_REF]. Furthermore, if f ′ ≥ 0, then this solution is unique.

We now turn to technical results, which will be used throughout the paper. The first one recalls some properties of Fokker-Planck type equations.

Lemma 2.2. Let b ∈ L ∞ (T d ; R d ), f ∈ L 2 (T d ), and g ∈ L 2 (T d ; R d ). There exists a weak solution m ∈ H 1 (T d ) to -∆m -div(mb) + λm = f + div(g) in T d .
Moreover, (i) m is the only element in L 2 (T d ) such that

(-∆ϕ + b • Dϕ + λϕ)m dx = f ϕ -g • Dϕ dx for every ϕ ∈ H 2 (T d ).
(ii) If f ≥ 0 and div(g) = 0, then, either f = 0 and m = 0, or m > 0.

(iii) If p ≥ 2, g ∈ L p (T d ; R d ), and f ∈ L q (T d ), with      q = dp/(d + p) if p = d and d ≥ 2, q > d/2 if p = d and d ≥ 2, q = 1 if d = 1, then m ∈ W 1,p (T d ) and there exists a positive constant C = C( b L ∞ , λ, d, p, q) such that m W 1,p ≤ C ( m L 1 + g L p + f L q ) . (6)
(iv) If g ∈ L p (T d ; R d ) and f ∈ L q (T d ) for p > d and q > d/2, then there exists a positive constant

C = C( b L ∞ , f L p/2 , g L p , λ, d, p, q) such that m L ∞ ≤ C. (7)
Proof. The existence of a weak solution m ∈ H 1 (T d ) (and its uniqueness in H 1 (T d )) is standard and follows from the Lax-Milgram theorem (see [START_REF] Evans | Partial Differential Equations[END_REF][START_REF] Gilbarg | Elliptic Partial Differential Equations of Second Order[END_REF]). For (i) let ξ ∈ L 2 (T d ) and let v ∈ H 2 (T d ) be the unique solution (again see [START_REF] Evans | Partial Differential Equations[END_REF][START_REF] Gilbarg | Elliptic Partial Differential Equations of Second Order[END_REF]) to

-∆v + b • Dv + λv = ξ in T d .
Then using m as a test function for the equation satisfied by v we have

ξm dx = (-∆v + b • Dv + λv)m dx = f v -g • Dv dx.
If f, g = 0, and since ξ is arbitrary, we obtain m = 0 which proves (i). When ξ ≥ 0, div(g) = 0 and f ≥ 0 we obtain m ≥ 0 since, in this case, v ≥ 0 by maximum principle (see [START_REF] Gilbarg | Elliptic Partial Differential Equations of Second Order[END_REF]Theorem 8.1]). The strict positivity when m = 0 is then a consequence of the Harnack inequality [START_REF] Gilbarg | Elliptic Partial Differential Equations of Second Order[END_REF]Theorem 8.20]. This proves (ii)

The W 1,p estimate (iii) is a direct consequence of [6, Theorem 1.7.4] while (iv) follows from De Giorgi-Nash-Moser estimates [START_REF] Gilbarg | Elliptic Partial Differential Equations of Second Order[END_REF]Theorem 8.17].

Remark 2.3. We recall from [START_REF] Leoni | A first course in Sobolev spaces[END_REF]Theorem 10.41] that, for 1 < r ≤ ∞, any element h ∈ W -1,r (T d ) can be identified with g 1 +div(g 2 ) for

g 1 ∈ L r (T d ), g 2 ∈ L r (T d ; R d ) and h W -1,r = ( g 1 r L r + g 2 r L r ) 1/r .
In particular the conclusions of Lemma 2.2 can be extended to equations with right-hand side in W -1,r (T d ) for appropriate values of r.

The following proposition contains a priori estimates on classical solutions to (1).

Proposition 2.4. Assume f ∈ C 0 b (R).
Then there exists a positive constant

K = K( f L ∞ , d) such that for every classical solution (u, m) to (1) it holds that u L ∞ ≤ f L ∞ λ , (8)
Du L ∞ ≤ K, (9)
and

∆u L ∞ ≤ 2 f L ∞ + K 2 2 =: M. ( 10
)
Furthermore, if λ > M , then m L ∞ ≤ λ λ -M m 0 L ∞ . (11)
Proof. Inequality (8) is a direct consequence of the comparison principle for u and ( 9) is given by [39, Theorem 1.1], since any classical solution is also a continuous viscosity solution. The estimate [START_REF] Bressan | Generic properties of first-order mean field games[END_REF] directly follows from the equation satisfied by u. For the last inequality [START_REF] Brezis | Functional analysis, Sobolev spaces and partial differential equations[END_REF], notice that we may rewrite the second equation in [START_REF] Achdou | Finite Difference Methods for Mean Field Games[END_REF] as

-∆m -Du • Dm + (λ -∆u)m = λm 0 in T d .
Therefore when λ > M we deduce from (10) that λ-∆u > 0 and therefore the Fokker-Planck equation satisfies a comparison principle from which the inequality follows.

Stable solutions

The following definition is taken from [START_REF] Briani | Stable Solutions in Potential Mean Field Game Systems[END_REF].

Definition 2.5 (Stable solutions). Let f ∈ C 1 b (R) and let (u, m) be a classical solution to [START_REF] Achdou | Finite Difference Methods for Mean Field Games[END_REF]. We say that (u, m) is stable if (v, ρ) = (0, 0) is the unique classical solution to

-∆v + Du • Dv + λv = f ′ (m)ρ in T d , -∆ρ -div (ρDu) + λρ = div (mDv) in T d . ( 12
)
Definition 2.6. Let f ∈ C 1 b (R) and let (u, m) be a classical solution to [START_REF] Achdou | Finite Difference Methods for Mean Field Games[END_REF]. A pair

(v, ρ) ∈ H 1 (T d ) × L 2 (T d ) is a weak solution to (12) if it satisfies Dv • Dϕ + Du • Dvϕ + λvϕ dx = f ′ (m)ρϕ dx for every ϕ ∈ H 1 (T d ) (13)
and

(-∆ψ + Du • Dψ + λψ)ρ dx = -mDv • Dψ dx for every ψ ∈ H 2 (T d ). ( 14
) Lemma 2.7. Let f ∈ C 1 b (R) ∩ C 1,1 loc (R) and (u, m) ∈ C 2,α (T d ) × C 2,α (T d ) be a classical solution to (1) and (v, ρ) ∈ H 1 (T d ) × L 2 (T d ) be a weak solution to (12). Then (v, ρ) ∈ C 2,α (T d ) × C 2,α (T d
) and is a classical solution to [START_REF] Brezzi | Finite Dimensional Approximation of Nonlinear Problems: Part I: Branches of Nonsingular Solutions[END_REF].

Proof. Notice first that if ρ ∈ C 0,α (T d ), then by Schauder estimates [START_REF] Gilbarg | Elliptic Partial Differential Equations of Second Order[END_REF]Corollary 6.3]

we have v ∈ C 2,α (T d ) and then also ρ ∈ C 2,α (T d ). It is therefore enough to prove that if (v, ρ) ∈ H 1 (T d ) × L 2 (T d ) is a weak solution to (12) then ρ ∈ C 0,α (T d ). For this, let 1 < p < ∞, depending on d and α, be such that Morrey's inequality yields W 1,p (T d ) ֒→ C 0,α (T d ). It is now enough to prove that v ∈ W 1,p (T d ) which implies ρ ∈ W 1,p (T d ) by Lemma 2.2 (iii).
Since ρ ∈ L 2 (T d ) and f ′ is bounded we deduce from standard elliptic regularity that v ∈ H 2 (T d ). Therefore in the case d = 1, 2 the fact that v ∈ W 1,p (T d ) directly follows from Sobolev's inequality. In the rest of the proof we assume d ≥ 3.

Since v ∈ H 2 (T d ) we have that div(mDv) = Dm • Dv + m∆v is an element of L 2 (T d ).
It is then well known that the second equation in [START_REF] Brezzi | Finite Dimensional Approximation of Nonlinear Problems: Part I: Branches of Nonsingular Solutions[END_REF] has a weak solution ρ ∈ H 2 (T d ). Using the property (i) in Lemma 2.2, we deduce that in fact ρ = ρ and therefore ρ ∈ H 2 (T d ). In the case where d = 3, 4 we have ρ ∈ L r (T d ) for every 2 ≤ r < ∞ and in particular ρ ∈ L p (T d ). Injecting this information in the equation satisfied by v we conclude that v ∈ W 2,p (T d ) (see [START_REF] Gilbarg | Elliptic Partial Differential Equations of Second Order[END_REF]Theorem 9.11]) and the conclusion follows. We may therefore assume that d ≥ 5 and in this case Sobolev's inequality yields ρ ∈ L q1 (T d ), where q 1 = 2d d-4 and, arguing as above, that v ∈ W 2,q1 (T d ). We now know that div(mDv) belongs to L q1 (T d ). One may keep this bootstrap argument going and conclude that either we can obtain that ρ ∈ W 2,d/2 (T d ) after a finite number of steps, or there exists a sequence of real numbers 2 ≤ q n < d/2, defined by

q n+1 = dqn d-2qn , q 0 = 2,
and such that ρ ∈ W 2,qn (T d ) for every n ≥ 1.

We claim that the latter case cannot happen. Indeed, if it were the case, we notice that qn+1 qn = d d-2qn > 1 so that the sequence is increasing. In particular q n ≥ 2 for every n and therefore qn+1 qn

> d d-4 . It follows that d 2 > q n ≥ 2 d d-4
n which yields a contradiction.

We must therefore have ρ ∈ W 2,d/2 (T d ) after a finite number of steps and hence, also, ρ ∈ L s (T d ) for every s ∈ [1, ∞). Using elliptic regularity one more time we have v ∈ W 1,p (T d ). This concludes the proof according to the discussion at the beginning of the argument.

The following proposition states that stable solutions to (1), although not unique in general, are isolated. The result is a straightforward adaptation of [13, Proposition 4.2] with weaker norms. We provide the proof in Appendix B.

Proposition 2.8 (Stable solutions are isolated). Let f ∈ C 1 b (R), let (u, m
) be a stable solution to (1), and let p = max{2, d}. Then there exists

R > 0 such that, if (ũ, m) = (u, m) is another classical solution to (1), then u -ũ W 1,p + m -m L 2 > R.
In the next theorem we establish the fact that for any λ > 0 outside of a finite set, every solution to the MFG system (1) is stable. Furthermore, we are able to quantify the upper bound on the values of λ for which this property may fail. This upper bound has the following form

Λ := max 2 m 0 L ∞ , K 2 2 + m 0 L ∞ f ′ L ∞ . ( 15
)
where the constant K is given in Proposition 2.4

Theorem 2.9.

Let f ∈ C 1 b (R).
Then there exists a finite, and possibly empty, subset Σ ⊂ (0, Λ], where Λ is given in [START_REF] Camilli | A Convergence Rate for the Newton's Method for Mean Field Games with Non-Separable Hamiltonians[END_REF], such that, if λ ∈ (0, +∞) \ Σ, then every classical solution to [START_REF] Achdou | Finite Difference Methods for Mean Field Games[END_REF]

is stable. Furthermore, if f ′ ≥ 0, then Σ = ∅.
Before proving Theorem 2.9 we establish three necessary technical lemmata.

Lemma 2.10. Assume that f ∈ C 1 b (R) and let g ∈ L 2 (T d ), h ∈ H -1 (T d ) and (u, m) be classical solution to (1). There exists λ > 0 such that if λ ≥ λ, σ ∈ [0, 1] and (v, ρ) ∈ H 1 (T d ) × L 2 (T d ) is a weak solution to -∆v + Du • Dv + λv = σ (f ′ (m)ρ + g) in T d , -∆ρ -div(ρDu) + λρ = div(mDv) + σh in T d , (16)
then there exists a positive constant C > 0 such that

v H 1 + ρ L 2 ≤ C ( g L 2 + h H -1 ) .
Proof. Notice first that from Lemma 2.2 and Remark 2.3 we in fact have ρ ∈ H 1 (T d ). Moreover, recall from Young's inequality that if a, b ≥ 0, then for every ǫ > 0 there exists a positive constant C ǫ such that ab ≤ ǫa 2 + C ǫ b 2 . Finally, we recall from Proposition 2.4 that Du L ∞ ≤ K and we can assume that λ is large enough so that

m L ∞ ≤ 2 m 0 L ∞ .
Using ρ as a test function for the second equation in ( 16) we have

|Dρ| 2 + λ |ρ| 2 dx = -ρDu • Dρ + mDv • Dρ dx + σ h, ρ H -1 ,H 1 ≤ Du L ∞ ρ L 2 Dρ L 2 + m L ∞ Dv L 2 Dρ L 2 + h H -1 ρ H 1 ≤ 2ǫ Dρ 2 L 2 + C ǫ Du 2 L ∞ ρ 2 L 2 + m 2 L ∞ Dv 2 L 2 + ǫ ρ 2 H 1 + C ǫ h 2 H -1 ≤ 3ǫ Dρ 2 L 2 + ǫ ρ 2 L 2 + C ǫ K 2 ρ 2 L 2 + 4 m 0 2 L ∞ Dv 2 L 2 + h 2 H -1 .
Therefore, if we choose ǫ = 1/3 and if λ > C 1/3 K 2 + 1/3 we obtain

ρ 2 L 2 ≤ C 1/3 λ -C 1/3 K 2 -1/3 -1 4 m 0 2 L ∞ Dv 2 L 2 + h 2 H -1 . (17)
We set

η := C 1/3 λ -C 1/3 K 2 -1/3 -1
. We now use v as a test function for the first equation in [START_REF] Cardaliaguet | Learning in Mean Field Games: The Fictitious Play[END_REF] to obtain the estimate

|Dv| 2 + λ |v| 2 dx = -vDu • Dv + σf ′ (m)ρv + σgv dx ≤ Du L ∞ Dv L 2 v L 2 + f ′ L ∞ ρ L 2 v L 2 + g L 2 v L 2 ≤ ǫ Dv 2 L 2 + ρ 2 L 2 + C ǫ Du 2 L ∞ v 2 L 2 + f ′ 2 L ∞ v 2 L 2 + 1 2 g 2 L 2 + v 2 L 2 ≤ ǫ 1 + 4η m 0 2 L ∞ Dv 2 L 2 + ǫη h 2 H -1 + 1 + C ǫ K 2 + f ′ 2 L ∞ v 2 L 2 + g 2 L 2 .
Hence, if we choose

ǫ < 1 + 4η m 0 2 L ∞ -1
and if λ is large enough we deduce that there exists a positive constant C such that v 2

H 1 ≤ C h 2 H -1 + g 2 L 2 .
The conclusion then follows from [START_REF] Cardaliaguet | The Master Equation and the Convergence Problem in Mean Field Games[END_REF].

Lemma 2.11. Assume that f ∈ C 1 b (R), let g ∈ L 2 (T d ), h ∈ H -1 (T d ),
and (u, m) be a classical solution to [START_REF] Achdou | Finite Difference Methods for Mean Field Games[END_REF]. There exists λ > 0 such that if λ ≥ λ there exists a weak solution

(v, ρ) ∈ H 2 (T d ) × H 1 (T d ) to -∆v + Du • Dv + λv = f ′ (m)ρ + g in T d , -∆ρ -div(ρDu) + λρ = div(mDv) + h in T d , (18)
and v H 2 + ρ H 1 ≤ C ( g L 2 + h H -1 ) .
Proof. The proof relies on the Leray-Schauder fixed point theorem [START_REF] Gilbarg | Elliptic Partial Differential Equations of Second Order[END_REF]Theorem 11.3]. Fix ρ ∈ L 2 (T d ).

Then there exists a unique solution v ∈ H 2 (T d ) to the equation

-∆v + Du • Dv + λv = f ′ (m)ρ + g in T d , and a unique ρ ∈ H 1 (T d ) satisfying -∆ρ -div(ρDu) + λρ = div(mDv) + h in T d
in the weak sense. In this way we define a mapping Φ :

L 2 (T d ) → L 2 (T d ) by setting Φ(ρ) = ρ.
This mapping is continuous and compact. In order to apply the Leray-Schauder theorem we have to prove that the solutions to ρ = σΦ(ρ) are bounded in L 2 (T d ) for every σ ∈ [0, 1]. This uniform bounded is provided by Lemma 2.10 since the equation ρ = σΦ(ρ) is equivalent to the system [START_REF] Cardaliaguet | Learning in Mean Field Games: The Fictitious Play[END_REF]. This proves that Φ has a fixed point which, by construction, provides a solution to [START_REF] Carmona | Probabilistic Theory of Mean Field Games with Applications I-II. Probability Theory and Stochastic Modelling[END_REF].

Lemma 2.12. Let f ∈ C 1 b (R) and (u, m) be a classical solution to [START_REF] Achdou | Finite Difference Methods for Mean Field Games[END_REF]. Then (u, m) is stable provided that one of the following conditions holds:

(i) (Monotonicity of the coupling) f ′ ≥ 0.
(ii) (Large discount factor) λ > Λ, where Λ is given in [START_REF] Camilli | A Convergence Rate for the Newton's Method for Mean Field Games with Non-Separable Hamiltonians[END_REF].

Proof. Let (v, ρ) be a weak solution to [START_REF] Brezzi | Finite Dimensional Approximation of Nonlinear Problems: Part I: Branches of Nonsingular Solutions[END_REF]. Using Lemma 2.7 we may assume that

(v, ρ) ∈ C 2 (T d ) × C 2 (T d ).
Using ρ as a test-function in [START_REF] Briani | Stable Solutions in Potential Mean Field Game Systems[END_REF] and v as a test-function in [START_REF] Caloz | Numerical analysis for nonlinear and bifurcation problems[END_REF] we get

Dv • Dρ + ρDu • Dv + λvρ dx = f ′ (m)ρ 2 dx, (19)
and

Dρ • Dv + ρDu • Dv + λρv dx = -m |Dv| 2 dx. (20)
Subtracting ( 20) from ( 19) we obtain

m |Dv| 2 dx = -f ′ (m) |ρ| 2 dx. (21)
In the case where f ′ ≥ 0, and since m is nonnegative (see Lemma 2.2), this implies that Dv = 0 on the set where m > 0, in particular mDv = 0 on T d . Using Lemma 2.2, we deduce that ρ = 0 on T d . Using the uniqueness of the solution to the equation satisfied by v, we conclude that v = 0. We now assume that λ > Λ. Notice that this and Proposition 2.4 yield m L ∞ ≤ 2 m 0 L ∞ . Then, using ρ as a test function in [START_REF] Caloz | Numerical analysis for nonlinear and bifurcation problems[END_REF], it follows from ( 9), the positivity of m and ( 21), that

|Dρ| 2 + λ |ρ| 2 dx ≤ Du L ∞ ρ L 2 Dρ L 2 + Dρ L 2 mDv L 2 ≤ K ρ L 2 Dρ L 2 + m 1/2 L ∞ Dρ L 2 m 1/2 Dv L 2 ≤ K ρ L 2 Dρ L 2 + m 1/2 L ∞ Dρ L 2 -f ′ (m) |ρ| 2 1/2 ≤ Dρ 2 L 2 + K 2 2 ρ 2 L 2 + 1 2 f ′ L ∞ m L ∞ ρ 2 L 2 ≤ Dρ 2 L 2 + K 2 2 + m 0 L ∞ f ′ L ∞ ρ 2 L 2 .
From our assumption on λ, we obtain ρ L 2 = 0 and the conclusion follows as in the first case.

Proof of Theorem 2.9. Notice first that the case λ > Λ as well as the case f ′ ≥ 0 directly follow from Lemma 2.12. Therefore, we only have to consider the case λ ∈ (0, Λ] with non-monotone coupling. Assume that

(v, ρ) ∈ H 1 (T d )× L 2 (T d
) is a nontrivial solution to [START_REF] Brezzi | Finite Dimensional Approximation of Nonlinear Problems: Part I: Branches of Nonsingular Solutions[END_REF]. We define the linear operator

B ∈ L(H 1 (T d ) × L 2 (T d ), L 2 (T d ) × H -1 (T d )) by B(w, µ) = (Du • Dw -f ′ (m)µ, -div(µDu + mDw)) .
Then, with obvious abuse of notations, (v, ρ) is a solution to [START_REF] Brezzi | Finite Dimensional Approximation of Nonlinear Problems: Part I: Branches of Nonsingular Solutions[END_REF] if and only if

(-∆ + λI + B) (v, ρ) = 0. ( 22
)
Let now λ be as in the conclusion of Lemma 2.10. We may assume that λ > Λ. Then ( 22) is equivalent to

-∆ + λI + B (v, ρ) = λ -λ (v, ρ). (23)
In other words, if (v, ρ) is a nontrivial solution to [START_REF] Brezzi | Finite Dimensional Approximation of Nonlinear Problems: Part I: Branches of Nonsingular Solutions[END_REF], then λ-λ belongs to ev(L), the set of eigenvalues of the linear operator

L := -∆ + λI + B, where L ∈ L(H 2 (T d ) × H 1 (T d ), L 2 (T d ) × H -1 (T d )). From
Lemmata 2.11 and 2.12, we know that L is invertible and, from the compact embedding of

H 2 (T d ) × H 1 (T d ) into L 2 (T d ) × H -1 (T d ), that L -1 is compact as a element of L(L 2 (T d ) × H -1 (T d ))
. It then follows from standard results on compact operators (see [START_REF] Brezis | Functional analysis, Sobolev spaces and partial differential equations[END_REF]Theorem 6.8] for instance) that the spectrum sp(L -1 ) of L -1 is at most countable and can only have an accumulation point at 0. Since, for µ = 0, µ ∈ ev(L) implies µ -1 ∈ sp(L -1 ), we deduce that the set ev(L) ∩ ( λ -Λ, λ) is finite. Hence there can only exist finitely many λ ∈ (0, Λ) such that ( 23) holds for a nontrivial solution (v, ρ) to [START_REF] Brezzi | Finite Dimensional Approximation of Nonlinear Problems: Part I: Branches of Nonsingular Solutions[END_REF].

Reformulation of the MFG system

We are now going to reformulate (1) as an abstract equation F (u, m) = 0. In order to do this let us introduce an additional assumption which will be useful to state some results in a general form.

(H) X and Z are Banach spaces such that

C 2,α (T d ) × C 2,α (T d ) ⊂ X ⊂ H 1 (T d ) × L 2 (T d ) and Z ⊂ L 1 (T d ) × (W 1,∞ (T d )) ′ . ( 24
)
with continuous embeddings. For every (v, ρ) ∈ X, we have

G(v, ρ) := 1 2 |Dv| 2 -f (ρ), -div (ρDu) -λm 0 belongs to Z (25)
and, for every (ξ, ζ) ∈ Z, the equation

-∆v + λv = ξ in T d , -∆ρ + λρ = ζ in T d , (26)
admits a unique distributional solution T (ξ, ζ) that belongs to X.

Assumption (H) allows to define two mappings G : X → Z and T : Z → X such that, at least formally, (u, m) solves (1) if and only

F (u, m) := (u, m) + T (G(u, m)) = 0. ( 27
)
This is made rigorous in the following result. ) if and only if it satisfies [START_REF] Fleming | Controlled Markov Processes and Viscosity Solutions[END_REF].

Proposition 3.1. Assume that f ∈ W 1,∞ (R) and that (H) holds. A pair (u, m) ∈ C 2,α (T d ) × C 2,α (T d ) solves (1
Proof. Assume that (u, m) ∈ C 2,α (T d ) × C 2,α (T d
) is a solution to [START_REF] Achdou | Finite Difference Methods for Mean Field Games[END_REF]. We can rewrite (1) as

-∆u + λu = f (m) -1 2 |Du| 2 in T d , -∆m + λm = λm 0 + div (mDu) in T d ,
or, with an obvious abuse of notation,

(-∆ + λI)(u, m) = -G(u, m). ( 28
)
We can apply T = (-∆ + λI) -1 on both sides of (28) and use the linearity of T to obtain that

(u, m) = -T (G(u, m)).
Conversely if (u, m) satisfies [START_REF] Fleming | Controlled Markov Processes and Viscosity Solutions[END_REF] we can apply the operator (-∆ + λI) to get [START_REF] Gilbarg | Elliptic Partial Differential Equations of Second Order[END_REF].

Let now (u, m) be a classical solution to (1) and X and Z be Banach spaces such that (H) holds and (u, m) ∈ X. Formally, the Fréchet differential of the mapping F defined in [START_REF] Fleming | Controlled Markov Processes and Viscosity Solutions[END_REF] at (u, m) is given by

dF [u, m](v, ρ) := (v, ρ) + T (dG[u, m](v, ρ)) (29)
where

dG[u, m](v, ρ) = (Du • Dv -f ′ (m)ρ, -div (ρDu) -div (mDv)) . (30)
The rigorous proof of Fréchet differentiability will be made in Section 4 below. Notice that (v, ρ) solves [START_REF] Brezzi | Finite Dimensional Approximation of Nonlinear Problems: Part I: Branches of Nonsingular Solutions[END_REF] if and only if

dF [u, m](v, ρ) = 0. ( 31
)
Theorem 3.2 (Isomorphism property of stable solutions). Let f ∈ C 1 b (R) and let X and Z be Banach spaces such that (H) holds and let Y be a Banach space such that

Y ⊂ X ⊂ H 1 (T d ) × L 2 (T d ) with Y ⊂ X compact.
Assume also that the mapping G : X → Z defined by ( 25) is continuously differentiable and that T ∈ L(Z, Y ). Then the mapping F : X → X defined by ( 27) is continuously differentiable with dF = I + T • dG, and, for every stable solution

(u, m) ∈ C 2,α (T d ) × C 2,α (T d ) to (1), the linear operator dF [u, m] is an isomorphism on X.
Proof. The differentiability of F is a direct consequence of the chain rule and the differentiability of

G. If (u, m) ∈ C 2,α (T d ) × C 2,α (T d ) is a stable solution to (1), then (I + T • dG[u, m]) is injective on C 2,α (T d ) × C 2,α (T d ). By Lemma 2.7 we know that any weak solution (v, ρ) ∈ H 1 (T d ) × L 2 (T d ) to (12) belongs to C 2,α (T d ) × C 2,α (T d ). Hence the operator (I + T • dG[u, m]) is also injective on H 1 (T d ) × L 2 (T d ), in particular it is injective on X. Notice now that since dG[u, m] ∈ L(X, Z) and T ∈ L(Z, Y ) we have that T • dG[u, m] ∈ L(X, Y ) and since Y ⊂ X is compact we deduce that the operator T • dG[u, m] ∈ L(X) is compact. Using
Fredholm's alternative, we conclude that dF [u, m] is an isomorphism on X.

We now provide two concrete examples of Banach spaces X, Y , and Z, satisfying (H) and the assumptions of Theorem 3.2. These examples will constitute the building blocks of the applications studied in Section 5.

Example 3.3. Let f ∈ C 1 b (R), let 0 < γ < β < α and set X = C 2,β (T d ) × C 2,γ (T d ), Y = C 2,α (T d ) × C 2,β (T d ) and Z = C 0,α (T d ) × C 0,β (T d )
where α is fixed in [START_REF] Bensoussan | Mean Field Games and Mean Field Type Control Theory[END_REF].

Let us check that (H) and the assumptions of Theorem 3.2 are satisfied in this case. From the Arzela-Ascoli theorem we have the compact embedding of Y into X. It is also clear that X ⊂ H 1 (T d )×L 2 (T d ).

We fix some stable solution

(u, m) ∈ C 2,α (T d ) × C 2,α (T d ) to (1). From the Schauder estimates [28, Corollary 6.3] we have T ∈ L(Z, Y ). It remains to prove that G : X → Z is well-defined. First notice that |Dv| 2 C 0,α ≤ 2 Dv 2 C 0,α ≤ C v 2 C 2,β .
Since f is Lipschitz continuous we also have

f (ρ) C 0,α ≤ f W 1,∞ ρ C 0,α ≤ C f W 1,∞ ρ C 2,γ
Finally, writting div(ρDv) = Dρ • Dv + ρ∆v, we obtain that

div(ρDv) C 0,α ≤ 2 Dρ C 0,α Dv C 0,α + 2 ρ C 0,α ∆v C 0,α ≤ C ρ C 2,γ v C 2,β .
In conclusion we have

G(v, ρ) Z ≤ C v 2 C 2,β + ρ 2 C 2,γ ≤ C (v, ρ) 2 X ,
and the mapping G is therefore well defined from X to Z.

In the case where G is also continuously differentiable we have dG : X → L(X, Z) and we can apply Theorem 3.

2 to conclude that dF [u, m] is an isomorphism on C 2,β (T d ) × C 2,γ (T d ). Example 3.4. Let f ∈ C 1 b (R) and X = W 1,p (T d ) × L q (T d ), Y = W 2,p/2 (T d ) × W 1,r (T d ) and Z = L p/2 (T d ) × W -1,r (T d )
where d < p, q < ∞ and r > 1 are such that 1/r ≥ 1/p + 1/q and we have a compact embedding W 1,r (T d ) ֒→ L q (T d ).

We now verify that (H) and the assumptions in Theorem 3.2 are satisfied. Since we assume p > d the Rellich-Kondrachov theorem gives the compact embedding W 2,p/2 (T d ) ֒→ W 1,p (T d ). It follows that there is a compact embedding from Y into X. Using W 2,p estimates [28, Theorem 9.11] and the W 1,p estimate (6) from Lemma 2.2 we have that

T ∈ L(Z, W 2,p/2 (T d ) × W 1,r (T d )). Therefore T ∈ L(Z, Y ) ⊂ L(Z, X). Moreover is is easy to check that G : X → Z is well defined.
Therefore in the case where G is continuously differentiable we obtain from Theorem 3.

2 that dF [u, m] is an isomorphism on W 1,p (T d ) × L q (T d ).
We now turn to the particular case of a monotone coupling and prove that the unique classical solution to [START_REF] Achdou | Finite Difference Methods for Mean Field Games[END_REF], which we also know to be stable by Theorem 2.9, satisfies a stronger isomorphism property than the one resulting from the direct application of Theorem 3.2.

Theorem 3.5 (Isomorphism property for monotone couplings). Let

f ∈ C 1 b (R) ∩ C 1,1 loc (R). Assume that f ′ ≥ 0 and let (u, m) ∈ C 2,α (T d ) × C 2,α (T d ),
where α is fixed in (5), be the unique classical solution to [START_REF] Achdou | Finite Difference Methods for Mean Field Games[END_REF]. Then (H) holds for X = C 2,α (T d )×C 2,α (T d ) and Z = C 0,α (T d )×C 0,α (T d ). Moreover, assume that G in (H) is continuously differentiable and define F according to [START_REF] Fleming | Controlled Markov Processes and Viscosity Solutions[END_REF]. Then dF [u, m] is an isomorphism on X.

Proof. According to Theorem 2.9, (u, m) is a stable solution to [START_REF] Achdou | Finite Difference Methods for Mean Field Games[END_REF]. Therefore the injectivity of dF [u, m] follows from the definition of stable solutions.

We now prove its surjectivity. More precisely we prove that for every (w, µ)

∈ C 2,α (T d )×C 2,α (T d ) there exists (v, ρ) ∈ C 2,α (T d ) × C 2,α (T d ) such that dF [u, m](v, ρ) = (w, µ), (32)
which is equivalent to showing that there exists a unique solution

(v, ρ) ∈ C 2,α (T d ) × C 2,α (T d ) to -∆v + Du • Dv + λv = f ′ (m)ρ -∆w + λw in T d , -∆ρ -div (ρDu) + λρ = div (mDv) -∆µ + λµ in T d . ( 33
)
The argument relies on the Leray-Schauder fixed point theorem and is adapted from [START_REF] Cardaliaguet | The Master Equation and the Convergence Problem in Mean Field Games[END_REF]Lemma 3.4]. Let us first define the mapping for which we will find a fixed point. Fix some ρ ∈ L 2 (T d ). Then there exists a solution v ∈ H 2 (T d ) to

-∆v + Du • Dv + λv = f ′ (m)ρ -∆w + λw in T d ,
and a unique weak solution ρ ∈ H 1 (T d ) to -∆ρ -div (ρDu) + λρ = div (mDv) -∆µ + λµ in T d . This allows us to define a mapping Φ : L 2 (T d ) → L 2 (T d ) by setting Φ(ρ) = ρ. In order to apply the Leray-Schauder fixed point theorem [START_REF] Gilbarg | Elliptic Partial Differential Equations of Second Order[END_REF]Theorem 11.3] we have to prove that the set of solutions ρ ∈ L 2 (T d ) to ρ = σΦ(ρ), where σ ∈ [0, 1], is uniformly bounded in L 2 (T d ). This amounts to prove a uniform bound on the solutions to

-∆v + Du • Dv + λv = σ [f ′ (m)ρ -∆w + λw] in T d , -∆ρ -div (ρDu) + λρ = div (mDv) -σ [∆µ -λµ] in T d .
Using ρ as a test function for the equation satisfied by v we obtain that

Dv • Dρ + ρDu • Dv + λvρ dx = σ f ′ (m)ρ 2 + ρ (λw -∆w) dx. (34)
Similarly, using v as a test function for the equation satisfied by ρ we obtain that

Dρ • Dv + ρDu • Dv + λρv dx = -m |Dv| 2 + σ (λvµ + Dv • Dµ) dx. (35)
Subtracting [START_REF] Lasry | Jeux à champ moyen. I -Le cas stationnaire[END_REF] to [START_REF] Huang | Large Population Stochastic Dynamic Games: Closed-Loop McKean-Vlasov Systems and the Nash Certainty Equivalence Principle[END_REF] we get

m |Dv| 2 dx = σ -f ′ (m)ρ 2 -ρ (λw -∆w) + λvµ + Dv • Dµ dx.
The positivity of f ′ then yields

m |Dv| 2 dx ≤ C ( ρ L 2 w H 2 + v H 1 µ H 1 ) . (36)
Let now ξ ∈ H -1 (T d ) and z ∈ H 1 (T d ) be the unique weak solution to

-∆z + Du • Dz + λz = ξ in T d .
Using z as a test function for the equation satisfied by ρ and recalling that m ≥ 0 (see Lemma 2.2-(ii)) we get

ξ, ρ H -1 ,H 1 = mDv • Dz + σ (λzµ + Dµ • Dz) dx ≤ C m 1/2 Dz L 2 m 1/2 Dv L 2 + z H 1 µ H 1 ≤ C z H 1 ( ρ L 2 w H 2 + v H 1 µ H 1 ) 1/2 + µ H 1 ≤ C ξ H -1 ( ρ L 2 w H 2 + v H 1 µ H 1 ) 1/2 + µ H 1
where we used [START_REF] Lasry | Jeux à champ moyen. II -Horizon fini et contrôle optimal[END_REF] to obtain the second inequality. Since ξ is arbitrary, we deduce by duality that

ρ H 1 ≤ C ( ρ L 2 w H 2 + v H 1 µ H 1 ) 1/2 + µ H 1 ,
and hence, by Young's inequality, we get

ρ H 1 ≤ C w H 2 + v 1/2 H 1 µ 1/2 H 1 + µ H 1 .
From elliptic regularity we know that

v H 2 ≤ C ( ρ L 2 + w H 2 ) ≤ C w H 2 + v 1/2 H 1 µ 1/2 H 1 + µ H 1 ,
and hence, after another application of Young's inequality, we get a uniform bound on v in

H 2 (T d ) v H 2 ≤ C ( w H 2 + µ H 1 ) .
Using standard H 1 estimates we can now deduce a uniform on ρ in H 1 (T d )

ρ H 1 ≤ C ( w H 2 + µ H 1 ) .
The compactness on Φ then follows from the Rellich-Kondrachov theorem and we can therefore apply the Leray-Schauder fixed point theorem to obtain a pair (v; ρ) ∈ H 2 (T d ) × H 1 (T d ) solving [START_REF] Brezzi | Finite Dimensional Approximation of Nonlinear Problems: Part I: Branches of Nonsingular Solutions[END_REF].

Using a bootstrap argument similar to the one used in Lemma 2.7, we can obtain that

(v, ρ) ∈ C 2,α (T d ) × C 2,α (T d ).

Differentiability of the mapping F

In this section we prove the differentiability of the mapping F defined by [START_REF] Fleming | Controlled Markov Processes and Viscosity Solutions[END_REF] in both Hölder and Sobolev spaces. This corresponds to the situations considered in Examples 3.3 and 3.4.

Differentiability in Hölder spaces

The following proposition summarizes the differentiability properties of the Nemytskii operator on Hölder spaces which we will use to prove the differentiability of the mapping G.

C 2,γ (T d ) ∋ u → H[u] ∈ C 0,β (T d ) is also continuously differentiable. Furthermore, if h ∈ C 2,1 loc (R), then dH is locally Lipschitz continuous in L(C 2,γ (T d ), C 0,β (T d )).
Proof. The first statement is proved in [START_REF] Nugari | Further Remarks on the Nemitskii Operator in Hölder Spaces[END_REF]Theorem 4.1] (this is where we need the C 2 assumption on f ). For the second one let us write J ∈ L(C 2,γ (T d ), C 0,β (T d )) for the natural injection of C 2,γ (T d ) into C 0,β (T d ). It is enough to notice that H • J is continuously differentiable from the chain rule.

For the local Lipschitz continuity let v ∈ C 2,γ (T d ). We have

(h ′ (u 1 ) -h ′ (u 2 ))v C 0,β ≤ v C 0,β h ′ (u 1 ) -h ′ (u 2 ) C 0,β ≤ v C 2,γ h ′ (u 1 ) -h ′ (u 2 ) C 0,β so that dH[u 1 ] -dH[u 2 ] L(C 2,γ ,C 0,β ) ≤ h ′ (u 1 ) -h ′ (u 2 ) C 0,β .
Therefore dH is locally Lipschitz continuous as soon as

C 2,γ (T d ) ∋ u → h ′ (u) ∈ C 0,β (T d )
is locally Lipschitzian, which is the case if h ′ ∈ C 1,1 loc according to [START_REF] Nugari | Further Remarks on the Nemitskii Operator in Hölder Spaces[END_REF]Theorem 3.1].

We are now going to prove the differentiability of the mapping F , defined by ( 27), in the situations described in Example 3.3 and Theorem 3.5.

Lemma 4.2. Let f ∈ C 1 b (R) ∩ C 2 (R) and 0 < γ ≤ β ≤ α, where α is set in (5). Set X = C 2,β (T d ) × C 2,γ (T d ) and Z = C 0,α (T d ) × C 0,β (T d ).
Then the mapping G : X → Z, defined by [START_REF] Ferreira | Existence of weak solutions to stationary mean-field games through variational inequalities[END_REF], is continuously differentiable and, for every (u, m) ∈ X, we have

dG[u, m](v, ρ) = (Du • Dv -f ′ (m)ρ, -div(ρDu) -div(mDv)) . Furthermore, if f ∈ C 2,1
loc (R), then dG is locally Lipschitz continuous.

Proof. From Proposition 4.1, we have that the differentials of

C 2,β (T d ) ∋ u → |Du| 2 2 ∈ C 0,β (T d ) and C 2,γ (T d ) ∋ m → f (m) ∈ C 0,β (T d ),
are given by

C 2,β (T d ) ∋ v → Du • Dv ∈ C 0,β (T d ) and C 2,γ (T d ) ∋ ρ → f ′ (m)ρ ∈ C 0,β (T d ),
respectively. The second component of dG being continuous and bilinear, its differentiability is also easy to check. The remaining conclusions follow from the assumptions on f and last assertion in Proposition 4.1.

Then fact that F is continuously differentiable when

f ∈ C 1 b (R) ∩ C 2 (R)
is then a direct consequence of Lemma 4.2 and the chain rule.

Proposition 4.3. Let f ∈ C 1 b (R) ∩ C 2 (R) and let 0 < γ ≤ β ≤ α, where α is fixed in (5). Set X = C 2,β (T d ) × C 2,γ (T d ) and Z = C 0,α (T d ) × C 0,β (T d ).
Then the mapping F : X → X defined by ( 27) is continuously differentiable with dF given by [START_REF] Girault | Finite Element Methods for Navier-Stokes Equations[END_REF]. Furthermore, if f ∈ C 2,1 loc (R), then dF is locally Lipschitz continuous in L(X).

Differentiability in Sobolev spaces

We begin with a preliminary result on the differentiability of the Nemytskii operator on Lebesgue spaces. Proposition 4.4. Let (Ω, A, µ) be a finite measure space and let 1 < p, q, r < ∞ with q < p and [START_REF] Lasry | Mean Field Games[END_REF] and define the Nemytskii operator H : L p (Ω, µ) → L q (Ω, µ) by

1/r = 1/q -1/p. Consider a function h ∈ C 1 (R) such that |h ′ (x)| ≤ C 1 + |x| p/r
H[u](x) = h(u(x)) for x ∈ Ω for µ-a.e.x ∈ Ω. Then H is continuously differentiable with dH[u](v) = h ′ (u)v. Proof. Let u, v ∈ L p (Ω, µ). Since h ∈ C 1 (R), we have h(u(x) + v(x)) -h(u(x)) = v(x) 1 0 h ′ (u(x) + tv(x)) dt
for µ-a.e. x ∈ Ω. Therefore, it follows from Holder's and Jensen's inequalities as well as Fubini's theorem that

H[u + v] -H[u] -h ′ (u)v q L q = v 1 0 h ′ (u + tv) -h ′ (u) dt q L q ≤ v q L p 1 0 h ′ (u + tv) -h ′ (u) dt q L r ≤ v q L p 1 0 h ′ (u + tv) -h ′ (u) q L r dt.
It follows that H is Fréchet differentiable at u if we are able to prove that

lim n→∞ sup t∈[0,1] h ′ (u + tv n ) -h ′ (u) L r = 0. ( 38
)
for every sequence v n in L p (Ω, µ) converging to 0. Let v n be such a sequence, we may assume that the convergence also holds in the µ-almost everywhere sense. Using the uniform continuity of h ′ on the set B(u(x), sup n |v n (x)|) we obtain that

lim n→∞ sup t∈[0,1] |h ′ (u(x) + tv n (x)) -h ′ (u(x))| r = 0 for µ-a.e. x ∈ Ω. (39)
For every σ > 1 we denote by c σ the smallest positive constant such that

(a + b) σ ≤ c σ (a σ + b σ ) for every a, b ≥ 0.
Using [START_REF] Lasry | Mean Field Games[END_REF] we have, for every t ∈ [0, 1],

|h ′ (u + tv n ) -h ′ (u)| r ≤ c r |h ′ (u + tv n )| r + |h ′ (u)| r ≤ c r C 1 + |u + tv n | p/r r + C 1 + |u| p/r r ≤ c 2 r C r (2 + |u + tv n | p + |u| p ) ≤ c 2 r C r (2 + c p |v n | p + (1 + c p ) |u| p ) so that sup t∈[0,1] |h ′ (u + tv n ) -h ′ (u)| r ≤ c 2 r C r (2 + c p |v n | p + (1 + c p ) |u| p ) µ-a.e. ( 40
)
Clearly the right-hand side in [START_REF] Lions | Cours Au College de France[END_REF] converges in L 1 (Ω; µ) and is therefore uniformly integrable according to the Lebesgue-Vitali theorem [START_REF] Bogachev | Measure Theory[END_REF]Theorem 4.5.4]. This proves that the left-hand side of ( 40) is uniformly integrable. Combining this last point with [START_REF] Ley | Lipschitz Regularity Results for Nonlinear Strictly Elliptic Equations and Applications[END_REF], we can use the other implication in the Lebesgue-Vitali theorem to establish lim

n→∞ sup t∈[0,1] |h ′ (u + tv n ) -h ′ (u)| r µ(dx) = 0. Since sup t∈[0,1] |h ′ (u + tv n ) -h ′ (u)| r µ(dx) ≤ sup t∈[0,1] |h ′ (u + tv n ) -h ′ (u)| r µ(dx)
we deduce [START_REF] Leoni | A first course in Sobolev spaces[END_REF]. This proves the differentiability of H.

For the continuity of dH, let u n be a sequence in L p (Ω, µ) converging to u ∈ L p (Ω, µ). Notice that

dH[u n ] -dH[u] L(L p ,L q ) ≤ h ′ (u n ) -h ′ (u) L r .
The continuity can then be proven using the Lebesgue-Vitali theorem similarly to what we did for differentiability.

We now come back to the differentiability of the mapping F for the situation described in Example 3.4.

Proposition 4.5. Let f ∈ C 1 b (R)∩C 1,1 loc (R), X = W 1,p (T d )×L q (T d ), and Z = L p/2 (T d )×W -1,r (T d )
, where d ≤ p, q < ∞, q > p/2, and 1/r ≥ 1/p + 1/q. Assume also that W 1,r (T d ) ֒→ L q (T d ). Then the mapping F : X → X defined by ( 27) is continuously differentiable.

Proof. Under these assumptions, we recall from Example 3.4 that (H) holds. From Proposition 4.4 we have that the mappings

L q (T d ) ∋ m → f (m) ∈ L p/2 (T d ) and W 1,p (T d ) ∋ u → |Du| 2 ∈ L p/2 (T d ) are continuously differentiable. Moreover W 1,p (T d ) × L q (T d ) ∋ (u, m) → div(mDu) ∈ W -1,r (T d )
is also continuously differentiable as a continuous bilinear operator. The differentiability of F then follows from the chain rule.

Applications

In this section we provide three applications of the isomorphism property of stable solutions to (1).

Stability under perturbations of the MFG system

In this section, our goal is to study perturbations of the mean field game system (1), or equivalently [START_REF] Fleming | Controlled Markov Processes and Viscosity Solutions[END_REF].

In the case of stable solutions, the isomorphism property obtained in Theorems 3.2 and 3.5 motivates the use of the implicit function theorem, allowing us to consider a large class of perturbations. In what follows we provide a simple example of this idea. Namely, we consider the system

-∆u + 1 2 |Du| 2 + λu = f (m) + ǫ f (m) in T d , -∆m -div(mDu) + λm = λ((1 -ǫ)m 0 + ǫm 1 ) in T d , (41)
where ǫ > 0 is a small parameter and f ∈ C 2 ((0, +∞)) and

m 1 ∈ C 0,α (T d ) ∩ P(T d ) are perturbations of f ∈ C 1 b (R) ∩ C 2 (R)
and m 0 , respectively. We recall that α is fixed in [START_REF] Bensoussan | Mean Field Games and Mean Field Type Control Theory[END_REF].

Let (u, m) ∈ C 2,α (T d ) × C 2,α (T d
) be a stable solution to (1) and set 0 < γ < β < α. Since m 0 = 0 we also have m > 0 from Lemma 2.2 (ii). In particular, there exists a bounded neighborhood O of m in C 2,γ (T d ) and η > 0 such that, for every m ∈ O, we have m ≥ η. Notice that if we define

E O := m(x) : m ∈ O, x ∈ T d , then we have that E O is bounded in R and inf E O ≥ η > 0. Since f ∈ C 2 ((0, +∞)) is Lipschitz continuous on E O , we deduce from Proposition 4.1 that the mapping O ∋ m → f ( m) ∈ C 0,β (T d ) is well-defined and continuously differentiable. Let X = C 2,β (T d ) × C 2,γ (T d ), Z = C 0,β (T d ) × C 0,γ (T d
) and consider the mappings G and T , defined by ( 25) and [START_REF] Fleming | Deterministic and Stochastic Optimal Control[END_REF], respectively. We recall that it was checked in Example 3.3 that (H) holds in this case. We introduce the mapping

Ĝ : C 2,β (T d ) × O × R + → C 0,β (T d ) × C 0,γ (T d ) defined by Ĝ(u, m, ǫ) = G(u, m) -ǫ( f (m), λ(m 1 -m 0 )).
Then, setting

F (u, m, ǫ) = Î + T • Ĝ (u, m, ǫ),
where T is defined by [START_REF] Fleming | Deterministic and Stochastic Optimal Control[END_REF] and Î(u, m, ǫ) = (u, m), it holds that (u ǫ , m ǫ ) ∈ X solves (41) if and only if

F (u ǫ , m ǫ , ǫ) = 0.
Arguing as in the proof of Proposition 4.3, we have that F is continuously differentiable on C 2,β (T d ) × O ×[0, +∞). Moreover, it follows from Theorem 3.2 and Example 3.

3 that d (u,m) F [u, m, 0] = dF [u, m]
is an isomorphism on X. We can apply the implicit function theorem to obtain the following result.

Proposition 5.1 (Sensitivity analysis).

Let f ∈ C 1 b (R) ∩ C 2 (R), f ∈ C 2 ((0, +∞)), m 1 ∈ C 0,α (T d ) ∩ P(T d ), 0 < γ < β < α, where α is fixed in (5), and (u, m) ∈ C 2,α (T d ) × C 2,α (T d ) be a stable solution to (1)
. Then, for some ǫ 0 > 0 and every ǫ ∈ [0, ǫ 0 ), there exists

(u ǫ , m ǫ ) ∈ C 2,β (T d ) × C 2,γ (T d ) such that F (u ǫ , m ǫ , ǫ) = 0, with (u ǫ , m ǫ ) = (u, m) -ǫdF (u, m) -1 T ( f (m), λ(m 1 -m 0 )) + o(ǫ), and
u ǫ -u C 2,β + m ǫ -m C 2,γ = O(ǫ).
Furthermore, (u ǫ , m ǫ ) is a stable solution to (41).

Remark 5.2. We may choose α = β = γ in Proposition 5.1 under the additional assumption that f ′ ≥ 0.

Finite Element approximation of the MFG system

Our goal here is to obtain error estimates for the finite element approximation of a stable classical solution to (1) by applying the following result of Brezzi-Rappaz-Raviart [START_REF] Brezzi | Finite Dimensional Approximation of Nonlinear Problems: Part I: Branches of Nonsingular Solutions[END_REF] (see also [START_REF] Girault | Finite Element Methods for Navier-Stokes Equations[END_REF]Section IV.3]).

Theorem 5.3 ([29, Theorem 3.3 and Remark 3.5]). Let V, W be Banach spaces, let T, T h ∈ L(W, V ), for every h > 0, and let G : V → W be a continuously differentiable mapping such that dG locally Lipschitz continuous in L(V, W ).

Set F = I + T • G and let x ∈ V be such that F (x) = 0 and dF [x] is an isomorphism on V . If lim h→0 T -T h L(W,V ) = 0, (42)
then there exists h 0 > 0 and a neighborhood O of x in V such that, for every 0 < h ≤ h 0 , there exists

x h ∈ V such that F h (x h ) := (I + T h • G)(x h ) = 0. (43)
Furthermore the following properties hold (i) dF h [x h ] is an isomorphism on V , (ii) x h ∈ O for every 0 < h ≤ h 0 and there is no other solution to [START_REF] Osborne | Analysis and Numerical Approximation of Stationary Second-Order Mean Field Game Partial Differential Inclusions[END_REF] in O, (iii) there exists a constant K > 0, independent of h, such that

x -x h V ≤ K (T -T h )G(x) V . (44)
We fix d ≤ 3 and consider the following Banach spaces

X = W 1,p (T d ) × L q (T d ) and Z = L p/2 (T d ) × H -1 (T d ), (45)
where p, q ∈ (3, 6) with 1 p + 1 q = 1 2 . Notice that under these assumptions we have q > p/2. Moreover the Rellich-Kondrachov theorem gives the compact embedding H 1 (T d ) ֒→ L q (T d ) since q < 6. Therefore all the assumptions in Example 3.4 are satisfied (with r = 2). In particular, the linear operator T ∈ L(Z, X), defined by [START_REF] Fleming | Deterministic and Stochastic Optimal Control[END_REF], and the mapping G : X → Z, given by [START_REF] Ferreira | Existence of weak solutions to stationary mean-field games through variational inequalities[END_REF], satisfy (H). We set r := dp p+d < min{p, d} and we notice that Sobolev's inequality implies that W 2,r (T d ) ֒→ W 1,p (T d ).

For every h > 0, let T h be a quasi-uniform family of periodic triangulations of [0, 1] d (see [START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF]Definition 4.4.13]). Let also V h ⊂ W 1,∞ (T d ) be the associated finite element space induced by P 1 Lagrange finite elements. We define

S h ∈ L(H -1 (T d ), H 1 (T d )) by S h ξ = v h where v h is the unique element in V h such that T d Dv h • Dφ h + λv h φ h dx = ξ, φ h H -1 ,H 1 for every φ h ∈ V h .
In addition, we denote by S ∈ L(H -1 (T d ), H 1 (T d )) the linear operator defined by Sξ = v, where v ∈ H 1 (T d ) is the unique weak solution to

-∆v + λv = ξ in T d . (46)
These linear operators are known to be well defined through the Lax-Milgram theorem and we have

S L(H -1 ,H 1 ) , S h L(H -1 ,H 1 ) ≤ 1 min{1, λ} . (47)
We also have from [START_REF] Ern | Theory and Practice of Finite Elements[END_REF]Theorem 3.16,Theorem 3.18] that

(S -S h )ξ L 2 ≤ Ch (S -S h )ξ H 1 ≤ Ch Sξ H 1 ≤ Ch ξ H -1 , (48)
where the constant C is independent of h. Let ξ ∈ L r (T d ). According to [START_REF] Ern | Theory and Practice of Finite Elements[END_REF]Theorem 3.21] and [START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF]Theorem 8.5.3], there exists h 0 > 0 and a positive constant C such that, for every h ≤ h 0 , it holds that

S h ξ W 1,p ≤ C Sξ W 1,p . (49)
Using the continuous embedding W 2,r (T d ) ֒→ W 1,p (T d ) and W 2,p estimates [28, Theorem 9.11], we have

Sξ W 1,p ≤ C Sξ W 2,r ≤ C ξ L r . (50)
Combining ( 49) and (50), for every h ≤ h 0 , we have

S h ∈ L(L r (T d ), W 1,p (T d ))
, where S h L(L r ,W 1,p ) is bounded from above by a constant which is independent of h. Moreover, up to the choice of a smaller h 0 , from [23, Theorem 3.21, Corollary 3.23] and [START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF]Theorem 8.5.3] we also have, for every h ≤ h 0 , that

S h ∈ L(L p (T d ), W 1,p (T d )) with (S -S h )ξ W 1,p ≤ Ch Sξ W 2,p ≤ Ch ξ L p , (51)
where the constant C is independent of h.

We can now define the linear operator

T h ∈ L(L r (T d ) × H -1 (T d ), W 1,p (T d ) × H 1 (T d )),
with range in X h := V h × V h , by setting

T h (ξ, ζ) = (S h ξ, S h ζ) for every (ξ, ζ) ∈ L r (T d ) × H -1 (T d ).
Notice that from ( 26), ( 48) and (49), we also have

T h ∈ L(L p (T d ) × H -1 (T d ), W 1,p (T d ) × L 2 (T d )) with (T -T h )(ξ, ζ) W 1,p ×L 2 ≤ Ch (ξ, ζ) L p ×H -1 . (52)
In particular, using Sobolev inequalities ([11, Theorem 9.9, Corollary 9.11, Theorem 9.12], we have obtained that

(T -T h ) ∈ L(L r (T d )×H -1 (T d ), W 1,p (T d )×L s (T d )))∩L(L p (T d )×H -1 (T d ), W 1,p (T d )×L 2 (T d )),
where s is the critical Sobolev exponent for the continuous embedding

H 1 (T d ) ֒→ L s (T d ), namely      s = 6 if d = 3, s ∈ [6, ∞) if d = 2, s = ∞ if d = 1.
We set θ = p-d p and θ ⋆ = (p-2)s-2p (s-2)p so that

2 p = 1 -θ r + θ p and 1 q = p -2 2p = 1 -θ ⋆ s + θ ⋆ 2 .
Using complex interpolation (see [START_REF] Lunardi | Interpolation Theory[END_REF]Chapter 2]) we have that

L p/2 (T d ) = L r (T d ), L p (T d ) θ and L q (T d ) = L s (T d ), L 2 (T d ) θ ⋆ .
It follows from [START_REF] Lunardi | Interpolation Theory[END_REF]Theorem 2.6] together with (49), (50), (51), ( 47) and ( 48) that that

S -S h L(L p/2 ,W 1,p ) ≤ S -S h 1-θ L(L r ,W 1,p ) S -S h θ L(L p ,W 1,p ) ≤ Ch θ and S -S h L(H -1 ,L q ) ≤ S -S h 1-θ ⋆ L(H -1 ,L s ) S -S h θ ⋆ L(H -1 ,L 2 ) ≤ Ch θ ⋆ . Noticing that T -T h = (S -S h , S -S h ), we deduce that T -T h L(L p/2 ×H -1 ,W 1,p ×L q ) ≤ S -S h L(L p/2 ,W 1,p ) + S -S h L(H -1 ,L q ) ≤ C h θ + h θ ⋆ . Since      θ ⋆ = (p -3)/p if d = 3, θ ⋆ ∈ [(p -3)/p, (p -2)/p) if d = 2, θ ⋆ = (p -2)/p if d = 1,
we have σ := min{θ, θ ⋆ } = θ ⋆ . Therefore we obtain that

T h , (T -T h ) ∈ L(Z, X)
with

T -T h L(Z,X) ≤ Ch σ , (53)
for every h ≤ h 0 ≤ 1.

We can now apply Theorem 5.3 to obtain the following result.

Theorem 5.4 (Local convergence of finite element approximations). Let f ∈ W 2,∞ (R) and let (u, m) be a stable solution to (1). Let X, Z be defined according to (45) with 3 < p, q < 6 and 1/2 = 1/p + 1/q and let T h be defined as above. There exists h 0 ∈ (0, 1] and a neighborhood O of the origin in X such that, for every 0 < h ≤ h 0 , there exists a solution (u h , m h ) ∈ X h to

F h (u h , m h ) := (u h , m h ) + T h (G(u h , m h )) = 0 satisfying (i) (u, m) -(u h , m h ) ∈ O, (ii) (u h , m h ) is the unique solution to F h (u h , m h ) = 0 in (u, m) + O, (iii)
There exists a positive constant C > 0 such that

(u -u h , m -m h ) X ≤ Ch σ , where      σ = (p -3)/p if d = 3, σ ∈ [(p -3)/p, (p -2)/p) if d = 2, σ = (p -2)/p if d = 1, (iv) dF h [u h , m h ] ∈ L(X) is an isomorphism.
Proof. Using Proposition 4.5 we have that F is continuously differentiable and from Example 3.4 we know that dF [u, m] is an isomorphism on X. Note that the assumption f ∈ W 2,∞ (R) ensures that dG is locally Lipschitz continuous in L(X, Z). Moreover, from (53) we deduce that ( 42) is satisfied. We can therefore apply Theorem 5.3 to obtain the conclusion.

As a direct consequence of Theorem 5.4 (iv) we deduce the local convergence of Newton's method for the discretized problem.

Corollary 5.5 (Local convergence of the discrete Newton method). Under the assumptions of Theorem 5.4, let (u h , m h ) ∈ X h be a solution to

F h (u h , m h ) = 0. ( 54
)
Then there exists a neighborhood O of (u h , m h ) such that, if (u 0 h , m 0 h ) ∈ O, then the sequence (u k h , m k h ) given by Newton's method applied to (54), i.e.,

(u k+1 h , m k+1 h ) + T h G(u k h , m k h ) + dG[u k h , m k h ](u k+1 h -u k h , m k+1 h -m k h ) = 0, (55)
converges quadratically to (u h , m h ) in X h . in T d , -∆ρ -div(ρDu k ) + λρ = div(m k Dv) -div(m k Du k ) in T d . (60)
In the case where f ′ ≥ 0, we may use Theorem 3.5 instead of Example 3.3 to obtain a slightly better result than Theorem 5.8.

Theorem 5.10. Let f ∈ C 1 b (R) ∩ C 2 (R) with f ′ ≥ 0. Let (u, m) ∈ C 2,α (T d ) × C 2,α (T d
) be a stable solution to [START_REF] Achdou | Finite Difference Methods for Mean Field Games[END_REF]. Then there exists a neighborhood O of of (u, m) in C 2,α (T d ) × C 2,α (T d ) such that if (u 0 , m 0 ) ∈ O then the sequence (u k , m k ) generated by Newton's method (59) converges super-linearly to (u, m) in C 2,α (T d ) × C 2,α (T d ). Furthermore, if we also assume f ∈ C 2,1 loc (R), then the convergence is quadratic.

Finally we may also set X = W 1,p (T d ) × L q (T d ) and Z = L p/2 (T d ) × W -1,r (T d ), where d < p, q < ∞ and r > 1 is such that 1/r = 1/p+1/q and large enough so that there is a compact embedding W 1,r (T d ) ֒→ L q (T d ). Then, we may replace Proposition 4.3 and Example 3.3 by Proposition 4.5 and Example 3.4, respectively, in the discussion above to obtain the convergence of Newton's method in Sobolev spaces. The point being that in this case the neighborhood for the initial guess is expected to be less restrictive.

Theorem 5.11. Let f ∈ C 1 b (R) ∩ C 1,1 loc (R)
, let X and Z be as above and (u, m) ∈ C 2,α (T d ) × C 2,α (T d ) be a stable solution to (1). Then there exists a neighborhood O of (u, m) in W 1,p (T d ) × L q (T d ) such that, if (u 0 , m 0 ) ∈ O, then the sequence (u k , m k ) generated by Newton's method (59) converges superlinearly to (u, m) in W 1,p (T d ) × L q (T d ). Furthermore, if f ∈ W 2,∞ (R), then the convergence is quadratic.

A Proof of Theorem 2.1

We are going to apply Schauder's fixed point theorem in C 0,β (T d ) for some β ∈ (0, α] to be determined.

Fix some m ∈ C 0,β (T d ). Since f is assumed to be Lipschitz continuous, we have that f (m) ∈ C 0,β (T d ). From the standard theory of elliptic equation (see [START_REF] Gilbarg | Elliptic Partial Differential Equations of Second Order[END_REF]Theorem 15.12], using the gradient bound in Proposition 2.4, for the result with Dirichlet boundary conditions) we know that there exists a unique classical solution u ∈ C 2,β (T d ) to

-∆v + 1 2 |Dv| 2 + λv = f (m) in T d . (61)
Then from standard Schauder theory (see [START_REF] Gilbarg | Elliptic Partial Differential Equations of Second Order[END_REF]Corollary 6.3]) and Lemma 2.2, we also have a unique solution m ∈ C 2,β (T d ) to

-∆ρ -div (ρDu) + λρ = λm 0 in T d . (62)
This defines a mapping Φ : C 0,β (T d ) → C 0,β (T d ) by setting Φ(m) = m. We now prove that Φ is continuous. Let (m n ) n≥0 be a sequence in C 0,β (T d ) converging to some m in C 0,β (T d ). In particular, this sequence is bounded in C 0,β (T d ). Using the fact that f is Lipschitz continuous it follows that f (m n ) is also bounded in C 0,β (T d ). From the inequality (9) in Proposition 2.4, we have the existence of a positive constant K such any solution classical solution u n to (61), with m replaced by m n , satisfies

Du n L ∞ ≤ K and the constant K depends on the right-hand side of (61) only through f L ∞ (and hence is independent of n). Then, from [28, Theorem 13.6] we deduce that there exist constants γ ∈ (0, 1) and K > 0, depending on K and independent of n, such that

Du n C 0,γ ≤ K.
Using Schauder estimates [28, Corollary 6.3], one has that the sequence u n of solutions to (61) associated to m n is bounded in C 2,β (T d ) for β = min{α, γ}. From the Arzela-Ascoli theorem, it admits a subsequence converging in C 2 (T d ) to a solution u to (61) associated to m. Since this solution is unique the whole sequence must converge to this limit u. Then, using again Schauder estimates we also have that Φ(m n ) is bounded in C 2,β (T d ) by a constant depending on K, and a similar argument shows that it must converge to Φ(m).

We now claim the Φ(C 0,β (T d )) is bounded in C 0,γ (T d ) for γ ∈ (β, 1). Indeed let us choose 1 < p < ∞, depending only on d and γ, such that W 1,p (T d ) ֒→ C 0,β (T d ). From the W 1,p estimates (6) in Lemma 2.2, we have that Φ(C 0,β (T d )) is bounded in W 1,p (T d ) and therefore also in C 0,γ (T d ). Using the Arzela-Ascoli theorem we deduce that Φ(C 0,β (T d )) is compact in C 0,β (T d ).

We can now apply Schauder's fixed point theorem [START_REF] Gilbarg | Elliptic Partial Differential Equations of Second Order[END_REF]Corollary 11.2] to obtain a classical solution to [START_REF] Achdou | Finite Difference Methods for Mean Field Games[END_REF].

The argument for uniqueness is a straightforward adaptation of the one introduced in [START_REF] Lasry | Mean Field Games[END_REF].

B Proof of Proposition 2.8

The

(u, m) in W 1,p (T d ) × L 2 (T d ).
Note that we may assume that the convergence also holds in the almost everywhere sense. We then set

δ n = (u n , m n ) -(u, m) W 1,p ×L 2 and v n = δ -1 n (u n -u), ρ n = δ -1
n (m n -m). Then for every n the pair (v n , ρ n ) is a classical solution to

-∆v n + λv n = g n in T d , -∆ρ n + λρ n = div(h n ) in T d , (63)
where

g n = δ -1 n f (m n ) -f (m) + 1 2 |Du| 2 - 1 2 |Du n | 2 and
h n = δ -1 n (m n Du n -mDu) . We then notice that

g n L 2 ≤ δ -1 n ( f ′ L ∞ m n -m L 2 + K Du n -Du L 2 ) ≤ f ′ L ∞ + K and h n L 2 ≤ δ -1 n (K m n -m L 2 + m L ∞ Du n -Du L 2 ) ≤ C
, where K is given in Proposition 2.4 and C is a constant independent of n. Applying standard elliptic regularity to (63) is follows that (v n , ρ n ) is bounded in H 2 (T d )×H 1 (T d ). From the Rellich-Kondrachov theorem we deduce that it converges, up to a subsequence, to some (v, ρ) in W 1,q (T d ) × L 2 (T d ) for every 2 ≤ q < 2d d-2 . We may also assume that (v n , ρ n ) converges weakly to (v, ρ) in H 2 (T d )× H 1 (T d ). Moreover v n is bounded in W 1,d (T d ) (recall that we consider the case p = d) so that we may also assume that it converges weakly to v in W 1,d (T d ). We now claim that the following fact holds We then decompose

       δ -1 n (f (m n ) -f (m)) -→ f ′ (m)ρ in L 2 (T d ), ( 2δ
1 0 f ′ (λm n + (1 -λ)m)ρ n -f ′ (m)ρ dλ L 2 ≤ ρ n 1 0 f ′ (λm n + (1 -λ)m) -f ′ (m) dλ L 2 + f ′ (m)(ρ n -ρ) L 2 . ( 66
)
Since we assume that f ′ is bounded, the convergence to zero of the second term in the right-hand side of the last inequality follows from the convergence of ρ n to ρ in L 2 (T d ). Recalling that we assume that d > 2, it follows from Sobolev's inequality that H 1 (T d ) embeds continuously in L 2d d-2 (T d ). Using Hölder's inequality we can write

ρ n 1 0 f ′ (λm n + (1 -λ)m) -f ′ (m) dλ L 2 ≤ C ρ n H 1 1 0 f ′ (λm n + (1 -λ)m) -f ′ (m) dλ L d . ( 67
)
Using Jensen's inequality and Fubini's theorem we have This, together with (65), (66), (67), and (68), proves the first line in (64). For the second line of (64), we write

(2δ n ) -1 |Du| 2 -|Du n | 2 -Du • Dv L 2 = 1 2 (Du + Du n ) • Dv n -2Du • Dv L 2 ≤ 1 2 Du • (Dv n -Dv) L 2 + 1 2 Du n • Dv n -Du • Dv L 2 ≤ 1 2 Du L ∞ Dv n -Dv L 2 + 1 2 Du n • (Dv n -Dv) L 2 + 1 2 (Du n -Du) • Dv L 2 ≤ K Dv n -Dv L 2 + 1 2 (Du n -Du) • Dv L 2 ,
where the constant K is given in Proposition 2.4. Using Hölder's inequality we have We can now pass to the limit in the weak formulation of (63) to obtain that (v, ρ) ∈ H 1 (T d )×L 2 (T d ) is a weak solution to [START_REF] Brezzi | Finite Dimensional Approximation of Nonlinear Problems: Part I: Branches of Nonsingular Solutions[END_REF]. Indeed, for ϕ ∈ H 1 (T d ) we have This proves that (v, ρ) is a weak solution to [START_REF] Brezzi | Finite Dimensional Approximation of Nonlinear Problems: Part I: Branches of Nonsingular Solutions[END_REF]. According to Lemma 2.7, (v, ρ) is a classical solution and, since (u, m) is assumed to be stable, we must have (v, ρ) = (0, 0). This contradicts the fact that (v, ρ) W 1,p ×L 2 = 1 and concludes the proof.

(Du n -Du) • Dv L 2 ≤ Du n -Du L d Dv L 2d d-2 . It follows that (2δ n ) -1 |Du| 2 -|Du n | 2 -Du • Dv

Proposition 4 . 1 .

 41 Assume that h ∈ C 2 (R) and let β ∈ (0, 1]. Then the Nemytskii operatorH : C 0,β (T d ) → C 0,β (T d ) defined by H[u](x) = h(u(x)) for x ∈ T d is continuously differentiable and dH[u](v) = h ′ (u)v.In particular, for every γ ∈ (0, 1] the mapping

  n ) -1 |Du| 2 -|Du n | 2 -→ Du • Dv in L 2 (T d ), δ -1 n (m n Du n -mDu) ⇀ ρDu + mDv weakly in L 2d d+2 (T d ; R d ).

1 0f

 1 the first line of (64). Since f ∈ C 1 b (R) we have that δ -1 n (f (m n ) -f (m)) = ρ n ′ (λm n + (1 -λ)m) dλ a.e. in T d .(65)

1 0fd ≤ 1 0f 1 0f

 111 ′ (λm n + (1 -λ)m) -f ′ (m) dλ d L ′ (λm n + (1 -λ)m) -f ′ (m) d L d dλ.(68)Using the fact that f ′ is bounded and continuous and that m n converges a.e. to m one may apply Lebesgue's convergence theorem to obtainlim n→∞ ′ (λm n + (1 -λ)m) -f ′ (m) d L d dλ = 0.

L 2 ≤d- 2 .

 22 K Dv n -Dv L 2 + Du n -Du L d Dv L 2d Since v ∈ H 2 (T d ) ֒→ W 1, 2dd-2 (T d ) the right-hand side converges to 0. For the third line of (64), we first notice thath n = δ -1 n (m n Du n -mDu) = m n Dv n + ρ n Du. Let ϕ ∈ L 2d d-2 (T d ; R d ) = L 2d d+2 (T d ; R d ) ′ . From convergence of m n to m in L 2 (T d) and Hölder's inequality it is easy to see that m n ϕ converges to mϕ inL d d-1 (T d ; R d ) = (L d (T d ; R d )) ′ . Since v n converges weakly to v in W 1,d (T d ), we have that Dv n converges weakly to Dv in L d (T d ; R d ). It follows that lim n→∞ ϕ, m n Dv n (L 2d d+2 ) ′ ,L 2d d+2 = lim n→∞ m n ϕ • Dv n dx = lim n→∞ m n ϕ, Dv n (L d ) ′ ,L d = mϕ, Dv (L d ) ′ ,L d = ϕ, mDv (L 2d d+2 ) ′ ,L 2d d+2 .Since ρ n converges to ρ in L 2d d-2 (T d ) we have that ρ n Du converges to ρDu in L 2d d-2 (T d ), and because 2d d-2 ≥ 2d d+2 , we also have the strong (and therefore also weak) convergence of ρ n Du to ρDu in L 2d d+2 (T d ). This proves the third line of (64) and concludes the proof of the claim.

  lim n→∞ Dv n • Dϕ + λv n ϕ dx = Dv • Dϕ + λvϕ dx and lim n→∞ g n ϕ dx = f ′ (m)ρϕ -Du • Dvϕ dx, so that Dv • Dϕ + Du • Dvϕ + λvϕ dx = f ′ (m)ρϕ dx for every ϕ ∈ H 1 (T d ).Let now ψ ∈ H 2 (T d ) and recall that this implies that Dψ ∈ L 2d d-2 (T d ; R d ). Then we have lim n→∞ (-∆ψ + λψ)ρ n dx = (-∆ψ + λψ)ρ dx and lim n→∞ h n Dψ dx = (ρDu + mDv) • Dψ dx, so that (-∆ψ + Du • Dψ + λψ)ρ dx = mDv • Dψ dx for every ψ ∈ H 2 (T d ).

  argument is adapted from[START_REF] Briani | Stable Solutions in Potential Mean Field Game Systems[END_REF] Proposition 4.2]. We only prove the case d > 2, the other ones being simpler, in particular in what follows we have p = d. Assume that the conclusion is false, then there exists a sequence (u n , m n ) of classical solutions to (1) converging to

Acknowledgement. This work was partially supported by the ANR (Agence Nationale de la Recherche) through the COSS project ANR-22-CE40-0010 and the Centre Henri Lebesgue ANR-11-LABX-0020-01. The third author was partially supported by KAUST through the subaward agreement ORA-2021-CRG10-4674.6.

Proof. Since X h ⊂ X is finite dimensional, and since dF h [u h , m h ] ∈ L(X h ) is injective on X h by Theorem 5.4-(iv), we have that dF h [u h , m h ] is also an isomorphism on X h . We conclude using standard results on Newton's method (see Theorem 5.7 below).

Remark 5.6. Relation (55) amounts to finding (v h , ρ h ) ∈ X h such that

for every (φ, ψ) ∈ X h .

Newton's method

We recall here the convergence results Newton's method, see [ Theorem 5.7 (Classical Newton method). Let X and Y be Banach spaces and F : X → Y be continuously differentiable. Let x ∈ X be such that F (x) = 0 and dF [x] ∈ L(X, Y ) is an isomorphism. Then there exists a neighborhood O of x in X such that the sequence defined by

converges superlinearly to x. Furthermore, if dF is locally Lipschitz continuous in L(X, Y ), then the convergence is quadratic.

Let 0 < γ < β < α, where α is given in [START_REF] Bensoussan | Mean Field Games and Mean Field Type Control Theory[END_REF], and set

so that (H) is satisfied and consider the mapping T, G, and F defined by ( 26), [START_REF] Ferreira | Existence of weak solutions to stationary mean-field games through variational inequalities[END_REF], and (27), respectively. From Proposition 4.3 we know that the mapping F is continuously differentiable with dF given by [START_REF] Girault | Finite Element Methods for Navier-Stokes Equations[END_REF]. Moreover, if we fix a stable solution (u, m) ∈ C 2,α (T d ) × C 2,α (T d ), then we know from Example 3.3 that dF [u, m] is an isomorphism on X.

A direct application of Theorem 5.7 yields the following theorem.

Theorem 5.8.

) be a stable solution to [START_REF] Achdou | Finite Difference Methods for Mean Field Games[END_REF]. Then there exists a neighborhood O of (u, m) in C 2,β (T d ) × C 2,γ (T d ) such that, if (u 0 , m 0 ) ∈ O, then the sequence (u k , m k ) generated by Newton's method applied to (27), i.e.,

Furthermore, if we also assume f ∈ C 2,1 loc (R), then the convergence is quadratic. Remark 5.9. At each iteration, the relation (59) amounts to solving the linear system