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

Mathieu Laurière
email: lauriere@princeton.edu

Pierre-Louis Lions

Optimal control of conditioned processes with feedback controls

We consider a class of closed loop stochastic optimal control problems in finite time horizon, in which the cost is an expectation conditional on the event that the process has not exited a given bounded domain. An important difficulty is that the probability of the event that conditionates the strategy decays as time grows. The optimality conditions consist of a system of partial differential equations, including a Hamilton-Jacobi-Bellman equation (backward w.r.t. time) and a (forward w.r.t. time) Fokker-Planck equation for the law of the conditioned process. The two equations are supplemented with Dirichlet conditions. Next, we discuss the asymptotic behavior as the time horizon tends to `8. This leads to a new kind of optimal control problem driven by an eigenvalue problem related to a continuity equation with Dirichlet conditions on the boundary. We prove existence for the latter. We also propose numerical methods and supplement the various theoretical aspects with numerical simulations.

Introduction

In optimal control, the goal is generally to control the dynamics of a system so as to minimize a certain cost or, equivalently, to maximize a certain profit. If the evolution of the system is stochastic, the cost is computed in expectation over the possible realizations of the randomness. In this work, we consider a stochastic optimal control problem in which the optimization is performed conditionally on the occurrence or the non occurrence of a certain event.

Hence, what follows can be seen as an attempt to model limited rationality: the agent determines her strategy but disregards the possibility of some event. Such a situation unfortunately occurs in politics, for example when the effects of some decisions on the global climate change are not taken into account. There may even be situations when the probability of the event that conditionates the strategy becomes smaller and smaller, i.e. in the example above, the leader keeps thinking that the situation is safe whereas the probability of disastrous events becomes larger and larger.

In what follows, we focus on the case when the event of interest is that the process stays inside a given bounded domain, but of course one may think about many other situations.

A quite important part of the material contained in the present work directly comes from the lectures given by the third author in Collège de France in November and December 2016, see [START_REF] Lions | Cours du Collège de France[END_REF]. However, we will also discuss theoretical aspects including some technical details that were not dealt with in [START_REF] Lions | Cours du Collège de France[END_REF] for lack of time, propose numerical methods and illustrate the main ideas by numerical simulations. In particular, the latter will shed some light on some open problems related to the long time behavior.

In the sequel, we consider a bounded domain Ω Ă R d with a smooth boundary denoted by BΩ. The closure of Ω in R d is denoted by Ω. For a time horizon T ą 0, Q T stands for the time-space cylinder p0, T q ˆΩ. Let 1 Ω be the indicator function of Ω, taking values in t0, 1u such that 1 Ω pxq " 1 if and only if x P Ω. Let W " pW t q tě0 be a standard d-dimensional Brownian motion. We consider feedback controls (also called closed loop or Markovian controls) in L 8 pp0, `8q ˆRd ; R d q, which are deterministic functions of t and x. For such a control b P L 8 pp0, `8q ˆRd ; R d q, let X b " pX b t q tě0 be a solution to the SDE

dX b t " ´bpt, X b t qdt `σdW t . (1)
We assume that the distribution of X b 0 is absolutely continuous with respect to Lebesgue measure with density p 0 , and that p 0 is a smooth nonnegative function supported in Ω. The first time X b exits Ω is denoted by τ b : τ b " inftt ą 0 : X b t R Ωu. To motivate the problem that will be studied below, let us first focus on a simple case. Let g : Ω Ñ R d be a continuous function and L : Ω ˆRd Ñ R be a strictly convex function, continuously differentiable in its second argument. The problem is to minimize over feedback controls b P L 8 pp0, `8q ˆRd ; R d q the total cost

ż T 0 E " LpX b t , bpt, X b t qq | τ b ą t ‰ dt `E " gpX b T q | τ b ą T ‰ (2)
where X b has been defined above.

We can rewrite the above problem as an optimal control problem driven by a Kolmogorov-Fokker-Planck (KFP) equation. The infinitesimal generator associated to X b is

L b ϕ " ´σ2 2 ∆ϕ `b ¨Dϕ. (3)
The distribution of X b t has a density m b ptq, and m b solves the KFP equation:

B t m b `Lb m b " 0, in p0, T q ˆRd ,
with initial condition m b p0q " p 0 on Ω, where

L b ϕ " ´σ2 2 ∆ϕ ´divpbϕq (4)
is the dual operator associated to L b .

Let us now introduce p b satisfying the Dirichlet problem

$ ' & ' % B t p b `Lb p b " 0, in Q T ,
p b " 0, on p0, T q ˆBΩ, p b p0, ¨q " p 0 , on Ω.

(

) 5
It can be checked that for any smooth f :

R d Ñ R, E " f pX b t q1 tτ b ątu ‰ " ż Ω f pxqp b pt, xqdx.
In particular, ş Ω p b pt, xqdx " Ppτ b ą tq represents the probability of "survival" until time t. To alleviate the notations, this quantity will sometimes be denoted simply by ş Ω p b ptq. Hence the running cost in (2) can be written as [START_REF] Agmon | Estimates near the boundary for solutions of elliptic partial differential equations satisfying general boundary conditions[END_REF].

E " LpX b t , bpt, X b t qq | τ b ą t ‰ " E " LpX b t ,
We stress that to obtain the latter formulation, it is important to consider only controls that are in feedback form. In the sequel we will consider a slightly more general class of problems, in which the costs are given by functionals acting on p b pt,¨q ş Ω pptq . Remark 1.1. From (5) it is clear that only the restriction of b to Q T matters in the latter optimal control problem.

We will see below that ş Ω p b ptq decays as t grows. This represents an important difficulty at least for two aspects:

• in numerical simulations, the division by very small quantities is a problem

• the asymptotic behavior of the problem as the horizon tends to `8.

The paper is organized as follows. In Section 2, we present the theory on finite horizon conditional control problems (as it was discussed in [START_REF] Lions | Cours du Collège de France[END_REF]), focusing on the case when the control is bounded in L 8 by a fixed constant M : the existence of a minimizer as well as the optimality conditions are discussed; the latter have the form of a forward-backward system coupling a Dirichlet problem involving a forward Fokker-Planck equation and a Dirichlet problem involving a Bellman equation, as in the theory of mean field games, see [START_REF] Lasry | Jeux à champ moyen. II. Horizon fini et contrôle optimal[END_REF][START_REF] Lasry | Jeux à champ moyen. I. Le cas stationnaire[END_REF][START_REF] Lasry | Mean field games[END_REF], and in the theory of mean field type control see [START_REF] Carmona | Mean field forward-backward stochastic differential equations[END_REF][START_REF] Carmona | Probabilistic theory of mean field games with applications[END_REF][START_REF] Bensoussan | Mean field games and mean field type control theory[END_REF]. We prove existence of solutions for this system. Next, we describe a possible asymptotic behavior when the time horizon becomes large: to the best of our knowledge, proving these asymptotics is a difficult open problem. Plugging the ansatz into the previously mentioned system of PDEs leads to an interesting new system, coupling a principal eigenvalue problem involving a non symmetric second order elliptic operator and a stationnary Belmann equation. In order to study the latter problem, we give useful facts on principal eigenvalue problems for non symmetric operators in Section 3. In Section 4, we prove existence for the system that has been introduced at the end of Section 2, by showing first that this system can be seen as the first order optimality conditions of an optimal control problem driven by a principal eigenvalue problem (this material was also discussed in [START_REF] Lions | Cours du Collège de France[END_REF]). Section 5 is independent from Sections 3 and 4: it deals with the passage to the limit M Ñ `8 (M is the parameter introduced in Section 2), and relies very much on stability results for weak solutions of Fokker-Planck and Hamilton-Jacobi equations. For the finite horizon problem, a finite difference method, reminiscent of that introduced in [START_REF] Achdou | Finite difference methods for mean field games[END_REF][START_REF] Achdou | Mean field games: numerical methods[END_REF] is proposed in Section 6 and numerical simulations are performed for one and two dimensional examples. In particular, the results are in agreement with the ansatz made at the end of Section 2. Finally, in Section 7, a numerical method is proposed for the stationary problem discussed in Section 4, and the results include the asymptotics when the horizon tends to infinity.

The finite horizon case 2.1 Existence of a minimizer and necessary optimality conditions

As above, let L : Ω ˆRd Ñ R be a continuous function, strictly convex and continuously differentiable in its second argument. Let Φ, Ψ : L 2 pΩq Ñ R be Fréchet differentiable functionals which are bounded from below. The gradients of Φ and Ψ at p P L 2 pΩq are respectively denoted by F rps and Grps, that is, F rps, Grps P L 2 pΩq and pF rps, qq L 2 pΩq " DΦrpspqq, pGrps, qq L 2 pΩq " DΨrpspqq for all q P L 2 pΩq. Let ě 0 be a fixed constant (we stress that unless otherwise specified, can be 0). In view of Remark 1.1, we introduce the following notations for the set of controls: B " L 8 pQ T ; R d q, and for every M ą 0, B M " tb P L 8 pQ T ; R d q : }b} L 8 ď M u is the subset of controls bounded by M . Then, for b P B, we consider the cost functional

Jpbq " ż T 0 ˆşΩ p b pt, xqLpx, bpt, xqqdx ş Ω p b ptq `Φ " p b pt, xq ş Ω p b ptq ˙d t `Ψ " p b pT, ¨q ş Ω p b pT q  ´ ln ˆżΩ p b pT q ˙, (6)
where p b solves [START_REF] Agmon | Estimates near the boundary for solutions of elliptic partial differential equations satisfying general boundary conditions[END_REF]. We are going to address the following optimal control problem:

minimize J on B M . ( 7
)
Remark 2.1. The last term in (6) has a different nature from the other ones because it does not depend on the conditional probability p b pT,¨q ş Ω p b pT q . If ą 0, it penalizes the decay of the probability of survival. This term will be helpful in Section 5.

We start with the following result. Remark 2.3. Due to the lack of convexity of J, we do not know if the minimizer is unique.

Proof. Let us first give some relevant information on [START_REF] Agmon | Estimates near the boundary for solutions of elliptic partial differential equations satisfying general boundary conditions[END_REF] and prove that Jpbq is well defined for all b P B.

Standard arguments ensure that for all b P B, there exists a unique weak solution p P L 2 p0, T ; H 1 0 pΩqq X W 1,2 p0, T ; H ´1pΩqq to [START_REF] Agmon | Estimates near the boundary for solutions of elliptic partial differential equations satisfying general boundary conditions[END_REF]. Maximum norms estimates for equations in divergence form (see [START_REF] Lieberman | Second order parabolic differential equations[END_REF]Corollary 9.10])) tell us that there exists a positive constant p which depends on p 0 , Ω, T and }b} L 8 such that

0 ď p ď p, a.e. in Q T . (8)
Therefore, if }b} 8 ď M , }p} L 8 is bounded by a constant which depends only on p 0 , Ω, T and M . From ( 8), ( 5) can be written Bp Bt ´∆p `div pBq " 0, in Q T ,

where B P L 8 pQ T ; R d q. Hölder estimates for the heat equation with a right hand side in divergence form (see [START_REF] Lieberman | Second order parabolic differential equations[END_REF]Theorem 6.33]) tell us that p P C α{2,α pQ T q, for some 0 ă α ă 1, and that }p} C α{2,α pQ T q is bounded by a constant which depends only on p 0 , Ω, T and }b} 8 . Moreover, using the nonnegativity of p, a parabolic version of Harnack inequality, see [START_REF] Lieberman | Second order parabolic differential equations[END_REF]Theorem 6.27], and the fact that ş Ω p 0 " 1, we can prove by contradiction that p ą 0 in p0, T s ˆΩ. This allows one to define Jpbq for all b P B.

Consider now a minimizing sequence pb n q n in B M and p n the solution of (5) corresponding to b n . From the above estimates, we may assume that, up to the extraction of a subsequence, b n á b in L 8 pQ T ; R d q weak * and that p n tends to p in CpQ T q. It is easy to prove that p is the unique solution to (5) corresponding to b. Using the assumptions on L, Φ and Ψ, we can also see from the latter convergence that Jpbq ď lim inf Jpb n q.

Hence, b achieves the minimum of J.

Remark 2.4. The map b Þ Ñ p, where p is the solution of (5), is locally Lipschitz continuous from L 8 pQ T ; R d q to L 2 p0, T ; H 1 0 pΩqq X W 1,2 p0, T ; H ´1pΩqq. Hence, b Þ Ñ Jpbq is locally Lipschitz continuous using the assumptions on L, Φ and Ψ.

Next, for M ą 0, we obtain first order necessary optimality conditions for [START_REF] Berestycki | The principal eigenvalue and maximum principle for second-order elliptic operators in general domains[END_REF]

ı pxq ş Ω pptq `c1 ptq, in Q T , u " 0, on p0, T q ˆBΩ, upT, xq " 1 ş Ω ppT q G " ppT, ¨q ş Ω ppT q  pxq `c2 pT q, in Ω, (11)
$ ' ' ' & ' ' ' % B t ppt, xq ´σ2 2 ∆ppt, xq ´div ˆppt, ¨qH ξ ˆ¨, ˆżΩ pptq ˙Dupt, ¨q˙˙p xq " 0, in Q T , p " 0, on p0, T q ˆBΩ, pp0, ¨q " p 0 , in Ω, (12)
and c 1 , c 2 : r0, T s Ñ R are defined by

c 1 ptq " ´żΩ ppt, xq ´L `x, H ξ px, `şΩ pptq ˘Dupt, xqq ˘`F " ppt,¨q ş Ω pptq ı pxq ¯dx `şΩ pptq ˘2 , c 2 pT q " ´ ş Ω ppT q ´şΩ ppT, xqG " ppT,¨q ş Ω ppT q ı pxqdx `şΩ ppT q ˘2 . (13)
Proof. Let b P B M and let p b be the corresponding solution of (5). Let δb P L 8 pQ T ; R d q be a small variation of b, and δp be the corresponding variation of p (δp depends on b and δb but we do not write explicitly this dependence to save notations). Neglecting higher order terms thanks to Remark 2.4,

$ ' & ' % B t δp ´∆δp ´divpδpbq " divpp b δbq, in Q T ,
δp " 0, on p0, T q ˆBΩ, δpp0, ¨q " 0, in Ω.

At leading order, the variation of the cost is , .

δJ " ´ ş Ω δppT, xqdx ş Ω ppT, xqdx `1 ş Ω p b pT q ż Ω G " p b pT, ¨q ş Ω p b pT q  pxqδppT, xqdx ´şΩ δppT q `şΩ p b pT q ˘2 ż Ω p b pT, yqG " p b pT, ¨q ş Ω p b pT q  pyqdy `ż T 0 $ & % ş Ω p b pt,
dt.

Let us consider the adjoint problem (for the same control b):

$ ' ' ' ' ' ' ' & ' ' ' ' ' ' ' % ´Bt u ´∆u `b ¨Du " Lp¨, bq `F " p b ş Ω p b ı ş Ω p b ´şΩ p b ´Lp¨, bq `F " p b ş Ω p b ış Ω p b ˘2 1 Ω , in Q T , u " 0, on p0, T q ˆBΩ, upT, xq " ´ 1 Ω pxq ş Ω p b pT q `G " p b pT,¨q ş Ω p b pT q ı pxq ş Ω p b pT q ´şΩ p b pT, yqG " p b pT,¨q ş Ω p b pT q ı pyqdy `şΩ p b pT q ˘2 1 Ω pxq, in Ω.
Let us denote by u b the solution to the above problem. Using the equation on u b , integration by parts, and the equation on δp, we see that, at leading order δJ "

ż T 0 " ş Ω p b L b p¨, bq ¨δb ş Ω p b `żΩ δp p´B t u b ´∆u b `b ¨Du b q * dt `żΩ δppT, xqu b pT, xqdx " ż T 0 " ş Ω p b L b p¨, bq ¨δb ş Ω p b `xB t δp ´∆δp ´divpδpbq, u b y * dt " ż T 0 " ş Ω p b L b p¨, bq ¨δb ş Ω p b ´żΩ p b δb ¨Du b * dt.
Hence a necessary condition for b " b opt to be a minimizer of J over B M is that

b opt pt, xq " H ξ ˆx, ˆżΩ p bopt ptq ˙Du bopt pt, xq ˙.
Plugging this control back into the PDEs for p bopt and u p b opt , we obtain the forward-backward KFP-HJB system (11)-( 12).

Remark 2.7. As a consequence of Theorems 2.2 and (2.6), we get the existence of weak solutions for system (11)-( 12) with p P C 0 pQ T q X L 2 p0, T ; H 1 0 pΩqq X W 1,2 p0, T ; H ´1pΩqq and u P C 0 pr0, T s; L 2 pΩqq X L 2 p0, T ; H 1 0 pΩqq. Remark 2.8. Setting, for µ P p0, `8q, q Hpx, µ, ξq " H px, µξq µ " max

bPR d ; |b|ďM ˆξ ¨b ´Lpx, bq µ ˙, (14)
the HJB equation (11) can be written as follows:

´Bt u ´∆u `q H ˆ¨, ż Ω pptq, Du ˙`ż Ω ppt, yq q H µ ˆy, ż Ω pptq, Dupt, yq ˙dy1 Ω " F " pptq ş Ω pptq ı ş Ω pptq ´żΩ ˆpptqF " pptq ş Ω pptq ş Ω pptq ˘2 1 Ω , (15)
because, if b opt is optimal for J over B M , then

q H µ ˆy, ż Ω pptq, Dupt, yq ˙" Lpy, b opt pt, yqq `şΩ pptq ˘2 1 Ω .
Remark 2.9. If, for example, Lpx, bq " f pxq `1 2 |b| 2 , then the maximizer in ( 14) is µξ if µ}ξ} ď M and M ξ }ξ} otherwise. Thus q Hpx, µ, ξq "

# 1 2 µ|ξ| 2 ´fpxq µ , if µ|ξ| ď M, M |ξ| ´fpxq µ ´M2
2µ , otherwise.

Heuristics about long time behavior

For T ą 0, let p pT q and u pT q denote respectively the solution of the FP equation ( 12) and the HJB equation [START_REF] Bogachev | Fokker-Planck-Kolmogorov equations[END_REF] on the time interval r0, T s, and let b pT q opt denote the associated optimal control. Numerical simulations provided in Section 7.2.2 below indicate that, at least in some cases, the following asymptotic behavior as T Ñ 8 may be expected: for some η P p0, 1{2q, for all t s.t. t{T P pη, 1 ´ηq,

p pT q pt, xq ş Ω p pT q pt, yqdy Ñ ppxq, ˆżΩ p pT q pt, yqdy ˙upT q pt, xq Ñ ũpxq, 1 ş Ω p pT q pt, yqdy d dt ˆżΩ p pT q pt, yqdy ˙Ñ ´λ, b pT q opt pt, xq Ñ bpxq.
Notice that, if we multiply respectively the HJB equation ( 11) by p and the Fokker-Planck equation ( 12) by u, we integrate in pt, T q ˆΩ, and we sum the resulting equations, then, using the definitions of c 1 ptq and c 2 pT q, all the terms but one cancel, and we obtain ż Ω upt, xqppt, xqdx " ´ , for all t P r0, T s.

Plugging the ansatz into the FP equation ( 12), we obtain:

$ ' ' ' ' & ' ' ' ' % ´σ2 2 ∆p ´div ppH ξ p¨, Dũqq " λp, in Ω, p " 0 on BΩ, p ě 0 in Ω, ż Ω p " 1. ( 16
)
This means that λ is the principal eigenvalue of the Dirichet problem associated with the stationary FP equation with drift ´b " ´Hξ p¨, Dũq.

Moreover, plugging the ansatz in the HJB equation [START_REF] Bogachev | Fokker-Planck-Kolmogorov equations[END_REF], we find:

$ ' ' ' ' & ' ' ' ' % ´σ2 2 ∆ũ `Hp¨, Dũq " λũ `F rps `c1 1 Ω , in Ω, ũ " 0 on BΩ, ż Ω ũpxqppxqdx " ´ , c 1 " ´`ş Ω p pL p¨, H ξ p¨, Dũqq `F rpsq dx ˘. (17)
We show in the sequel that ( 16)-( 17) are the necessary optimality conditions of an optimization problem driven by an eigenvalue problem. This will yield existence for ( 16)- [START_REF] Evans | Partial differential equations[END_REF].

In [START_REF] Lions | Cours du Collège de France[END_REF], the third author discussed some mathematical tools that may be used in order to prove that the ansatz given above actually describes the asymptotic behavior as T Ñ 8, such as suitable generalizations to the notion of ergodic problems and to principal eigenvalue problems. However, we stress the fact that, to the best of our knowledge, proving rigorously the ansatz is an open problem. The main reason for that is that other asymptotic behaviors, involving for example time-periodic solutions, may be possible.

3 Facts about a principal eigenvalue problem

The principal eigenvalue related to an elliptic equation: known facts

We recall that Ω is a bounded domain of R d with a smooth boundary. Let b be a vector field: b P L 8 pΩ; R d q.

Recall that the elliptic operator L b is defined by [START_REF] Achdou | Mean field type control with congestion II : An augmented Lagrangian method[END_REF]. The weak maximum principle holds for L b , see [START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF]Theorem 8.1]: If v P H 1 pΩq is a weak subsolution of the Dirichlet problem

L b v " 0 in Ω v " 0 on BΩ i.e. is such that ż Ω Dv ¨Dw `wb ¨Dv ď 0,
for all w P H 1 0 pΩq such that w ě 0 almost everywhere in Ω, and if v `P H 1 0 pΩq, then v ď 0 almost everywhere in Ω.

It is also well known, see [START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF]Theorem 8.3], that for any g P L 2 pΩq, there exists a unique weak solution of the Dirichlet problem:

L b u " g in Ω, (18)
u " 0 on BΩ,

i.e. u P H 1 0 pΩq such that for all w P H 1 0 pΩq, ż

Ω Du ¨Dw `wb ¨Du " ż Ω gw. Moreover, u P H 2 pΩq X H 1 0 pΩq, see for example [17, §6.3, Theorem 4] and the operator g Þ Ñ u is bounded from L 2 pΩq to H 2 pΩq X H 1 0 pΩq. The operator S : L 2 pΩq Ñ L 2 pΩq defined by Spgq " u is compact. Let L 2
`pΩq be the closed cone made of the nonnegative functions in L 2 pΩq: from the weak maximum principle, we know that if g P L 2

`pΩq, then Sg P L 2 `pΩq.

From Agmon-Douglis-Nirenberg theorems, see [START_REF] Agmon | Estimates near the boundary for solutions of elliptic partial differential equations satisfying general boundary conditions[END_REF] or [START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF]Theorem 9.15], we know that if g P L p pΩq, for 2 ă p ă 8, then u P W 2,p pΩq X W 1,p 0 pΩq and the operator g Þ Ñ u is bounded from L p pΩq to W 2,p pΩq X W 1,p 0 pΩq. Let C 0 pΩq be the space of real valued continuous functions defined on Ω and vanishing on BΩ. For any 0 ă γ ă 1 that we fix, it will be useful to define the Banach space F " C 1,γ pΩq X C 0 pΩq and the cone F `" tv P F : v ě 0 P Ωu. Using the latter regularity results, we know that if g P CpΩq, then u P W 2,p pΩqXW 1,p 0 pΩq for all p ă 8. Sobolev embeddings of W 2,p pΩq X W 1,p 0 pΩq in F for p large enough imply that the linear operator g Þ Ñ u is compact from CpΩq to F. Let T be the restriction of the latter operator to F. Finally, as a consequence of Hopf lemma, which may be applied since b P L 8 pΩ; R d q, we know that if 0 ı g ě 0 in Ω, then Tg P IntpF `q. The previously mentioned facts allowes us to state the following theorem: Theorem 3.1. There exists a positive real number λ b ą 0, such that 1. 1{λ b is an eigenvalue of T 2. The related eigenspace is one-dimensional and can be written Ru where u P F `, u ą 0 in Ω 3. For all complex number µ such that µTpvq " v, for some non identically zero function v whose real and imaginary part belong to F, Repµq ě λ b .

Proof. When b P C 8 p Ω; R d q, Theorem 3.1 is exactly [17, Theorem 3, page 340]. However, in view of the facts that we have recalled above, the proof of [17, Theorem 3, page 340], which only requires the compactness of T, the strong maximum principle and Hopf lemma, can be repeated in the case when b P L 8 pΩ; R d q.

Remark 3.2. Theorem 3.1 implies that the spectral radius of T is positive and that it coincides with 1{λ b .

The positive number λ b is called the principal eigenvalue related to the Dirichlet problems ( 18)- [START_REF] Krein | Linear operators leaving invariant a cone in a Banach space[END_REF]. Setting e b " u{ ş Ω udx, we see that e b P F `is the unique weak solution of

´∆e b `b ¨De b " λ b e b in Ω, ( 20
)
e b " 0 on BΩ, (21) ż
Ω e b " 1, (22)
and e b ą 0 in Ω. If µ " 0 is an eigenvalue of S and w is any related eigenvector, then the regularity results for the Dirichlet problem

L b w " 1 µ w in Ω,
w " 0 on BΩ, imply that w P F, see for example [START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF]Theorem 9.15]. Hence, the spectrum and the eigenspaces of S and T coincide.

In [START_REF] Berestycki | The principal eigenvalue and maximum principle for second-order elliptic operators in general domains[END_REF], Berestycki, Nirenberg and Varadhan propose a definition of the principal eigenvalue of L b which can be applied even if BΩ is not smooth:

λ b " suptλ : Dφ P W 2,d loc pΩq, φ ą 0 in Ω, such that L b φ ´λφ ě 0 in Ωu. (23)
As stated in [START_REF] Berestycki | The principal eigenvalue and maximum principle for second-order elliptic operators in general domains[END_REF] (see also [START_REF] Nussbaum | On variational principles for the generalized principal eigenvalue of second order elliptic operators and some applications[END_REF]), if BΩ is smooth, then the definition in ( 23) is equivalent to the former definition of the principal eigenvalue of L b . Hence, if BΩ is smooth, ( 23) is a characterization of the principal eigenvalue of L b . Note that no restriction is made in [START_REF] Lehoucq | ARPACK users' guide, volume 6 of Software, Environments, and Tools[END_REF] on the behavior of φ near BΩ.

The main interest of ( 23) is that it is a variational characterization of the principal eigenvalue (recall that the standard characterization for self-adjoint operators is irrelevant in the present context): indeed, as proved in [START_REF] Berestycki | The principal eigenvalue and maximum principle for second-order elliptic operators in general domains[END_REF],

λ b " sup φPW 2,d loc pΩq inf Ω L b φ φ . ( 24
)
Note that in [START_REF] Donsker | On the principal eigenvalue of second-order elliptic differential operators[END_REF], Donsker and Varadhan gave a slightly different variational formula for the principal value of L b and proved [START_REF] Lieberman | Second order parabolic differential equations[END_REF] in the case when b is continuous and BΩ is smooth.

As a consequence of [START_REF] Lehoucq | ARPACK users' guide, volume 6 of Software, Environments, and Tools[END_REF], it is proven in [START_REF] Berestycki | The principal eigenvalue and maximum principle for second-order elliptic operators in general domains[END_REF] that for any domain Ω 1 Ă Ω, λ b pΩq ď λ b pΩ 1 q and that the inequality is strict if the set inclusion is strict (see Theorem 3.5 below, which is in fact a much stronger result).

An important consequence of ( 24), of the latter observation and of the maximum principle is an upper bound on λ b : Lemma 3.3 (Lemma 1.1. in [START_REF] Berestycki | The principal eigenvalue and maximum principle for second-order elliptic operators in general domains[END_REF]). Suppose that Ω contains a ball of radius R ď 1. Then

0 ă λ b ď C R 2
where C depends only on d and }b} L 8 pΩ;R d q .

The dependence of λ b with respect to b is also addressed in [START_REF] Berestycki | The principal eigenvalue and maximum principle for second-order elliptic operators in general domains[END_REF]:

Proposition 3.4 (Proposition 5.1 in [7]). The map b Þ Ñ λ b is locally Lipschitz continuous from L 8 pΩ; R d q to R.
The Lipschitz constant of the restriction of the latter map to B M " tb : }b} L 8 pΩ;R d q ď M u can be chosen to depend only on Ω and M . where δ depends only on Ω, }b} L 8 pΩ;R d q and δ 0 .

Some spectral properties of the adjoint operator

From Fredholm alternative, see [START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF]Theorem 5.11], we know that S and S ˚have the same spectra, and that if µ P C, dim pkerpI ´µSqq " dim pkerpI ´μS ˚qq.

Therefore 1{λ b is the spectral radius of S ˚, and λ b can be seen as the principal eigenvalue related to the dual boundary value problem with a Fokker-Planck equation: there exists a unique p b P H 1 0 pΩq such that

L b p b " λ b p b in Ω, ( 25
)
p b " 0 on BΩ, (26) ż
Ω p b " 1, ( 27
)
where L b is the adjoint elliptic operator defined by (4). Moreover, from Krein-Rutman theorem (see [START_REF] Krein | Linear operators leaving invariant a cone in a Banach space[END_REF] or e.g. [16, Theorem 1] applied to S, the eigenvector p b belongs to L 2 `pΩq.

Proposition 3.6. The normalized eigenvector p b characterized by ( 25)-( 27) is positive in Ω and continuous in Ω.

Proof. For b P L 8 pΩ; R d q, we now consider the Dirichlet problem related to the adjoint elliptic operator

L b L b q " g in Ω, (28)
q " 0 on BΩ.

We say that L b (respectively L b ) satisfies property (M) if for every v P H 1 0 pΩq, the fact that xL b v, wy is ě 0 (respectively xL b v, wy ě 0) for all nonnegative functions w P H 1 0 pΩq, implies that v is nonnegative in Ω. Property (M) for L b is a particular version of the weak maximum principle, so we have seen in paragraph 3.1 that it is true. On the other hand, since Ω has a smooth boundary, using Theorem 2.2.1. in [11, page 60], we see that property (M) is true for L b if and only if it is true for L b . Hence, L b enjoys property (M). Then, from Proposition 2.1.4. in [11, page 57], we know that for any g P H ´1pΩq, there exists a unique weak solution of the Dirichlet problem (28)-(29), i.e. a unique function q P H 1 0 pΩq such that for all w P H 1 0 pΩq, ż Ω Dq ¨Dw `qb ¨Dw " xg, wy.

Let us recall Hölder regularity results for the Dirichlet problem (28)-(29). First, using Theorem 8.15 in [START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF], we know that for s ą d{2, there exists a constant C which depends on d, Ω, }b} 8 and s such that if g P L s pΩq, then q P L 8 pΩq and }q} L 8 pΩq ď C}g} L s pΩq .

Next, from Theorem 8.29 in [START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF] and the latter L 8 bound, we see that if g P L s pΩq for s ą d{2, then there exists α P p0, 1q such that p P C α p Ωq and

}q} C α p Ωq ď C}g} L s pΩq ,
for a given constant C.

Let us now focus on the case when for s ą d{2, the function g belongs to L s pΩq and is nonnegative. The solution q of (28)-( 29) is nonnegative in Ω by property (M), and continuous. It is a weak supersolution of L b q " 0 in Ω. Hence, we can apply the weak Harnack inequality stated in Theorem 8.18 in [START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF]: for 1 ď t ă d{pd ´2q, for all y P Ω and R ą 0 such that B 4R pyq Ă Ω,

R ´d{t }q} L t pB 2R pyqq ď C inf B R pyq q,
for a constant C which depends on d, }b} 8 , R, s, and t. This implies that if q vanishes in Ω, then q must be identically 0 in Ω. Hence, if g is non identically 0, q must be positive in Ω.

Fix one of the Hölder exponents α obtained above and define the Banach space G " C α 0 pΩq and the cone G `" tg P G; g ě 0 P Ωu. The linear operator U : g Þ Ñ q is compact from G to G. Krein-Rutman theorem (the first version stated in paragraph 3.1) can be applied to U, which yields that there exists an eigenvector of U in G `associated with the spectral radius of U. Since the eigenvectors of U are eigenvectors of S, the latter must coincide with p b . Hence p b P Cp Ωq.

Finally, we can apply Harnack inequality to p b : p b is positive in Ω. Proof. We consider p n the eigenvector of L n related to λ n and normalized in L 2 pΩq norm: it solves

A stability result

´∆p n ´divpp n b n q " λ n p n in Ω, (30)
p n " 0 on BΩ, (31) ż
Ω p 2 n " 1. ( 32
)
We know that p n ą 0 in Ω. We know that p bn " p n { ş Ω p n . From Lemma 3.3, we know that 0 ă λ n is bounded from above by a constant independent of n. From this observation, it is easy to see that p n is bounded in H 1 0 pΩq by a constant independent of n. Extracting a subsequence, we may assume that p n tends to q weakly in H 1 0 pΩq, strongly in L 2 pΩq and that λ n tends to µ in R. We may pass to the limit and we get that ´∆q ´divpqbq " µq in Ω, q " 0 on BΩ, ż Ω q 2 " 1,

q ě 0, in Ω.
We also see that ş Ω p n tends to ş Ω q. The latter cannot be zero, otherwise q would be identically 0 in contradiction with ş Ω q 2 " 1. Hence ş Ω p n is bounded from below by a positive constant. Therefore p bn " p n { ş p n Ñ q " q{ ş Ω q weakly in H 1 0 pΩq and strongly in L 2 pΩq. If µ " λ b , then, using e b as a test function in the latter boundary value problem, we get ş Ω qe b " 0 which is not possible since 0 ı q ě 0 and e b ą 0. Therefore µ " λ b and q " p b . Since the two cluster points are unique, the whole sequences pλ n q and pp bn q tend to λ b and p b respectively.

Dependence of p b with respect to b

Lemma 3.8. The map b Þ Ñ p b is locally Lipschitz continuous from L 8 pΩ; R d q to H 1 0 pΩq.

Proof. First, using a similar argument as in the proof of Lemma 3.7, we see that the map b Þ Ñ p b is locally bounded from L 8 pΩ; R d q to H 1 0 pΩq. Take M ą 0 and two vector fields b 1 and b 2 such that }b 1 } L 8 pΩ;R d q ď M and }b 2 } L 8 pΩ;R d q ď M . Let us set δb " b 1 ´b2 , δp " p b1 ´pb2 , δλ " λ b1 ´λb2 :

´∆δp ´divpb 1 δpq " λ b1 δp `divpp b2 δbq `δλp b2 in Ω, (33)
δp " 0 on BΩ,

Ω δp " 0. ( (34) ż
) 35
Testing ( 33) by e b1 , we see that

xdivpp b2 δbq `δλp b2 , e b1 y H ´1pΩq,H 1 0 pΩq " ´żΩ p b2 δb ¨De b1 dx `δλ ż Ω p b2 e b1 dx " 0. ( 36
)
On the other hand, from Fredholm's theory, we know that for all g P H ´1pΩq such that xg, e b1 y H ´1pΩq,H 1 0 pΩq " 0, there exists z P H 1 0 pΩq a weak solution of

´∆z ´divpb 1 zq ´λb1 z " g in Ω, (37)
and that z is unique up to the addition of a multiple of p b1 . Moreover, the operator A : g Þ Ñ z is bounded from tg P H ´1pΩq, xg, e b1 y H ´1pΩq,H 1 0 pΩq " 0u to H 1 0 pΩq{pRp b1 q. Therefore, for all g P H ´1pΩq such that xg, e b1 y H ´1pΩq,H 1 0 pΩq " 0, there exists a weak solution z P H 1 0 pΩq of (37) such that }z} H 1 0 pΩq ď 2}A} }g} H ´1 pΩq . Then q " z ´pş Ω zqp b1 P H 1 0 pΩq is the unique weak solution of (37) such that ş Ω q " 0. We deduce that }q} H 1 0 pΩq ď C}g} H ´1pΩq for a positive constant C (independent of g). In what follows, the positive constant C (independent of g) will vary from one line to the other. From (33)-( 35), (36) and the previous arguments, we deduce that

}δp} H 1 0 pΩq ď C}divpp b2 δbq `δλp b2 } H ´1 pΩq ď C `}p b2 δb} L 2 pΩ;R d q `|δλ|}p b2 } L 2 pΩq ď C `}δb} L 8 pΩ;R d q `|δλ| ď C}δb} L 8 pΩ;R d q .
The third line is obtained by using the bound on }p b2 } L 2 pΩq discussed in the very first lines of the proof. The last line is obtained by using Proposition 3.4.

4 Long time behavior: an optimal control problem driven by a principal eigenvalue problem

and λ b is the principal eigenvalue related to (20)-( 22), (40) p b is the unique solution of ( 25)-( 27).

(41)

It will be convenient to use the following notations for the set of controls: r B " L 8 pR d ; R d q, and for every M ą 0, r B M " t b P L 8 pR d ; R d q : } b} L 8 ď M u is the subset of controls uniformly bounded by M .

Existence of an optimal solution

Lemma 4.1. For every M ą 0, there exists a minimizer of (38) subject to (40)-(41) over r B M .

Proof. We consider a minimizing sequence p bm q mą0 , } bm } 8 ď M and L m v " ´∆v `b m ¨Dv. Let λ bm be the principal eigenvalue of L m (with Dirichlet conditions) and let p bm be the unique solution of ( 25)-( 27) where b is replaced by bm .

We can extract a subsequence, still indexed by m, such that bm Ñ b in L 8 pΩ; R d q weakly *. Then from Lemma 3.7, λ bm tends to λ b ą 0 and p bm tends to p b in L 2 pΩq strongly and in H 1 0 pΩq weakly. Moreover since L is convex and continuous in its second argument and since bm á b in L 8 pΩ; R d q weak * and p bm converges to p b in L 2 pΩq, we know that

Necessary optimality conditions

We now derive a necessary optimality condition for the above problem, as we did in Theorem 2.6 in the finite time horizon problem. Theorem 4.2. Let M ą 0. If bopt P r B M is a minimizer of r J over r B M , then necessarily, for every x P R d , bopt pxq achieves the maximum in [START_REF] Boccardo | Existence de solutions non bornées pour certaines équations quasi-linéaires[END_REF] with ξ replaced by Dupxq where pp, uq solves the PDE system (16)- [START_REF] Evans | Partial differential equations[END_REF]. The necessary optimality condition can be written bopt pxq " H ξ px, Dupxqq .

Proof. Assume that b " bopt is a minimizer of (38) subject to (40)-(41). Let δ b P L 8 pΩ; R d q be a small variation of b, and δλ b, δp b be the corresponding variations of λ b and p b. We have already seen that

´∆δp b ´divpδp bb q " δλ bp b `λb δp b `divpp bδ bq in Ω, ( 42
)
δp b " 0 on BΩ, (43) ż
Ω δp b " 0. ( 44
for all q P H 1 0 pΩq. We see that the right hand side of (46) vanishes for q " p b. From Fredholm's theory, this implies the existence of a solution u and its uniqueness up to the addition of a multiple of e b.

Hence, at leading order,

For this choice of u b, using (42)-(44) yields

δ r J " ż Ω δ bpxqp bpxq ¨´L bpx, bpxqq ´Du bpxq ¯dx, (48)
and since p b ą 0 in Ω from Proposition 3.6, the first order optimality condition is that almost everywhere, δ bpxq ¨´L bpx, bpxqq ´Du bpxq ¯for all admissible variation δ b of b.

Let us recall that the Hamiltonian H is defined by [START_REF] Boccardo | Existence de solutions non bornées pour certaines équations quasi-linéaires[END_REF]. The first order optimality condition is equivalent to saying that almost everywhere in Ω, bpxq achieves the maximum in Hpx, Du bpxqq " max

ηPR d ; |η|ďM `Du bpxq ¨η ´Lpx, ηq ˘,
i.e., assuming that H is differentiable with respect to its second argument, bpxq " H ξ px, Du bpxqq.

(49)

We thus obtain the PDE system ( 16)- [START_REF] Evans | Partial differential equations[END_REF].

5 Back to the finite horizon problem, letting M Ñ 8

For simplicity only, we restrict ourselves to the following minimization problem

minimize J on B M , ( 50
)
where

Jpbq " ż T 0 ˜şΩ p b pt, xq `1 2 |bpt, xq| 2 `f pxq ˘dx ş Ω p b ptq ¸dt ´ ln ˆżΩ p b pT q ˙, (51)
and p b solves [START_REF] Agmon | Estimates near the boundary for solutions of elliptic partial differential equations satisfying general boundary conditions[END_REF]. Here, f is a smooth nonnegative function defined in Ω and the distribution p 0 of X 0 satisfies p 0 lnpp 0 q P L 1 pΩq. We also suppose that ą 0: this assumption will ensure a lower bound for ż Ω ppt, xqdx ą 0, @t P r0, T s, (58)

ż Q T |bpt, xq| 2 dxdt `żQ T ppt, xq|bpt, xq| 2 dxdt ă 8, (59)
ż Ω ppt, xqupt, xq " ´ for almost all t P p0, T q (60)

2. up to the extraction of a sequence, (a) p M Ñ p in L 1 pQ T q, almost everywhere, in W 1,q p0, T ; W ´1,q pΩqq weakly for all 1 ď q ď d`2 d`1 , (b) u M Ñ u in L 1 pQ T q, almost everywhere in Q T , and in L 2 p0, T ; H 1 0 pΩqq weakly, 3. (a) p is a distributional solution of ( 5), i.e for all ψ P C 8 pQ T q, such that ψ " 0 on r0, T s ˆBΩ and at t " T ,

ż Q T p b ˆBt ψpt, xq ´σ2 2 ∆ψpt, xq `bpt, xq ¨Dψpt, xq ˙dxdt " ż Ω p 0 pxqψp0, xqdx (61)
(b) u is a distributional solution of $ ' ' ' & ' ' ' % ´Bt upt, xq ´σ2 2 ∆upt, xq `H ``ş Ω pptq ˘Dupt, xq şΩ pptq " f pxq ş Ω pptq ´γptq, in Q T , upT, xq " ´ 1 Ω pxq ş Ω ppT, yqdy ,
in Ω (62)

4.

Jpbq ď lim

M Ñ`8 inf b 1 PB M Jpb 1 q " inf M ą0 inf b 1 PB M Jpb 1 q. ( 63
)
Remark 5.3. Equation (57) plays the role of the energy identity that is usually found in the theory of mean field games, see for example [START_REF] Porretta | Weak solutions to Fokker-Planck equations and mean field games[END_REF].

Besides, it seems possible to prove that γptq "

ż Ω ppt, xq `1 2 |bpt, xq| 2 `f pxq ˘dx `şΩ pptq ˘2 ,
by using the crossed regularity results in the spirit of what did A. Porretta for weak solutions of mean field games, but this would be too long for the present paper.

Proof. We start by looking for estimates uniform in M and compactness results.

For all what follows, it will be important to have a bound from below on ş Ω p M ptq uniform in M ; first, since f is bounded from below, and since ´ ln `şΩ p M pT q ˘ď Jpb M q ď Jp0q, we see that ş Ω p M pT q is bounded from below by a positive constant independent of M . Second, from known results on weak solutions of Fokker-Planck equations, [27, Proposition 3.10], we know that t Þ Ñ ş Ω p M ptq is a nonincreasing function from r0, T s to p0, `8q. Therefore ş Ω p M ptq is bounded from below by a positive constant uniform in t and M (but not in ), and from above by ş Ω p 0 .

Using this information and again the fact that Jpb M q ď Jp0q, we see that }γ M } L 1 p0,T q is bounded uniformly with respect to M . From the latter fact, we immediately deduce a bound on }u M } L 8 pQ T q uniform in M , by constructing subsolutions and supersolutions of the form pt, xq Þ Ñ θptq. Then testing (53) by u M e λu 2 M with λ large enough, we deduce that }u M } L 2 p0,T ;H 1 0 pΩqq is bounded uniformly with respect to M . This implies that ´Bt u M ´∆u M is bounded in L 1 pQ T q uniformly with respect to M . From classical results on the heat equation with L 1 data, see e.g. [START_REF] Boccardo | Nonlinear elliptic and parabolic equations involving measure data[END_REF][START_REF] Boccardo | Nonlinear parabolic equations with measure data[END_REF], we deduce that pu M q M is relatively compact in L 1 pQ T q and that there exists u P L 2 p0, T ; H 1 0 pΩqq X L 8 pQ T q such that, up to the extraction of a sequence, u M Ñ u in L 1 pQ T q and almost everywhere in Q T , even in L q pQ T q for all 1 ď q ă `8 and in L 8 pQ T q weak * due to the L 8 pQ T q bound, and Du M Ñ Du almost everywhere in Q T . Moreover t Þ Ñ upt, ¨q has a trace at t " T which is a Radon measure. Let us go back to the Fokker-Planck equation. First, the uniform bounds on }γ M } L 1 p0,T q and on ş

Ω p M ptq imply that ş Q T p M |b M | 2
is uniformly bounded with respect to M . Second, we deduce from (55) and the bound on }u M } L 2 p0,T ;H 1 0 pΩqq that }b M } L 2 pQ T q is bounded uniformly in M . The latter two points enable us to use [START_REF] Porretta | Weak solutions to Fokker-Planck equations and mean field games[END_REF]Proposition 3.10], which ensures that for all 1 ď q ă d`2 d`1 and 1 ď r ă d`2 d , there exists a constant C ą 0 uniform in M such that }p M } L 8 p0,T ;L 1 pΩqq `}p M } L r pQ T q `}Dp M } L q pQ T q `}B t p M } L q p0,T ;W ´1,q pΩqq ď C, where C depends on the bound on ş Q T p M |b M | 2 and on }p 0 } L 1 pΩq . Moreover, since p 0 P LlogLpΩq then for 1 ď q ď d`2 d`1 and 1 ď r

ď d`2 d , }p M } L 8 p0,T ;LlogLpΩqq `}p M } L r pQ T q `}Dp M } L q pQ T q `}B t p M } L q p0,T ;W ´1,q pΩqq ď C,
where C depends on the bound on [START_REF] Porretta | Weak solutions to Fokker-Planck equations and mean field games[END_REF]Proposition 3.10] ensures also that p M lies in a relatively compact subset of L 1 pQ T q. Therefore, there exists a nonnegative function p P L 1 pQ T q X W 1,q p0, T ; W ´1,q pΩqq for all 1 ď q ď d`2 d`1 , such that up to an extraction of a subsequence, we may assume that p M Ñ p in L 1 pQ T q, almost everywhere, and in W 1,q p0, T ; W ´1,q pΩqq weakly. Therefore, up to a further extraction, t Þ Ñ ş Ω p M pt, yqdy tends to t Þ Ñ ş Ω ppt, yqdy strongly in L 1 p0, T q and a.e. in p0, T q. Moreover, we know that p M belongs to Cpr0, T s; L 1 pΩqq and is bounded in L 8 p0, T ; LlogLpΩqq by a constant C independent of M . Let us fix t P r0, T s. There exists a sequence t n which tends to t such that }p M pt n q} LlogLpΩq ď C and p M pt n , ¨q Ñ p M pt, ¨q in L 1 pΩq. Up to the extraction of a subsequence, we can assume that p M pt n , xq Ñ p M pt, xq almost everywhere in Ω. Then Fatou lemma implies that ş Ω p M pt, xq logpp M pt, xqqdx ď lim inf tnÑt ş Ω p M pt n , xq logpp M pt n , xqqdx ď C. This implies that p M ptq is bounded in LlogLpΩq uniformly with respect to M and t. Hence, for all t, px Þ Ñ p M pt, xqq M is equiintegrable. By Dunford-Pettis theorem, p M pt, ¨q is relatively weakly compact in L 1 pΩq. Now choosing t " T , since we already know that p M pT, ¨q á ppT, ¨q weakly in W ´1,q pΩq, we deduce that p M pT, ¨q á ppT, ¨q weakly in L 1 pΩq.

ş Q T p M |b M | 2 and }p 0 } LlogLpΩq . Since }b M } L 2 pQ T q is bounded uniformly in M , the last point in
Let us now consider the sequence b M : from (55) and the almost everywhere convergence of Du M to Du and of p M to p, we deduce that b M tends to b almost everywhere in Q T , where bpt, xq " ˆżΩ pptq ˙Dupt, xq.

From Fatou's lemma and the uniform estimate on }b M } L 2 pQ T q , we see that b P L 2 pQ T q. Similarly, we see that ş Q T p|b| 2 is finite. We now use the uniform bound on

ş Q T p M |b M | 2 and the strong convergence of p M to p in L 1 pQ T q: for all measurable subset E of Q T , ż E p M |b M | ď ˆżQ T p M |b M | 2 ˙1 2 ˆżE p M ˙1 2 , which implies that p M b M is equiintegrable. Therefore p M b M á pb in L 1 pQ T q.
Hence, p is a distributional solution of (5), i.e., for all test function φ P C 8 pQ T q such that φ " 0 on r0, T sˆBΩ and on tT u ˆΩ, ż

Q T p´B t φ ´∆φ `b ¨Dφqp ´żΩ p 0 pxqφp0, xqdx " 0.
Then using [27, Theorem 3.6], p P Cpr0, T s; L 1 pΩqq and p is a renormalized solution of (5), because b P L 2 pQ T q and ş Q T p|b| 2 is finite. For φ as in Remark 5.1, we can pass to the limit in (56), because u M á u in L 8 pQ T q weak * and p M Ñ p in L 1 pQ T q, and we get

ż T 0 φ 1 ptq ż Ω upt, xqppt, xqdx " 0, which implies that t Þ Ñ ş Ω upt, xqppt, xqdx is constant. From (57), we know that ş Q T p M pt, xqu M pt, xqdxdt "
´ T . Passing to the limit, we see that

ş Q T ppt, xqupt, xqdxdt " ´ T . Since ş Ω upt, xqppt, xqdx does not depend on t, we conclude that ż Ω upt, xqppt, xqdx " ´ . (64)
We are left with passing to the limit in the Bellman equation ( 53). Since 0 ď γ M is bounded in L 1 p0, T q, ş T t γ M psqds is bounded in BV, and there exists a bounded positive measure γ on r0, T s such that up to a further extraction of a subsequence, γ M á γ in the sense of measures and µ M : µ M ptq " ş T t γ M psqds converges to some µ in L q p0, T q for all 1 ď q ă 8 and weakly * in L 8 p0, T q, and dµ dt " ´γ in the sense of distributions. We deduce that µ has bounded variations, so µpT q " lim tÑT µptq exists. For a smooth function ψ compactly supported in p0, T s, we get by passing to the limit that

ż T 0 ψptqdγptq ´ż T 0 ψ 1 ptqµptqdt " 0. (65)
Let us set v M pt, xq " u M pt, xq `µM ptq: we see that v M tends to v: vpt, xq " upt, xq `µptq in L q pQ T q for all 1 ď q ă 8, weakly * in L 8 p0, T q, and in L 2 p0, T ; H 1 pΩqq weakly. The function v M is a solution of the following boundary value problem:

$ ' ' ' ' ' ' & ' ' ' ' ' ' % ´Bt v M pt, xq ´σ2 2 ∆v M pt, xq `HM ``ş Ω p M ptq ˘Dv M pt, xq şΩ p M ptq " f pxq ş Ω p M ptq , in Q T , v M ´µM " 0, on p0, T q ˆBΩ, v M pT, xq " ´ 1 Ω pxq ş Ω p M pT, yqdy , in Ω. (66)
By using stability results for weak solutions of Bellman equations, see [START_REF] Boccardo | Existence de solutions non bornées pour certaines équations quasi-linéaires[END_REF][START_REF] Achdou | Convergence of a finite difference scheme to weak solutions of the system of partial differential equations arising in mean field games[END_REF], we obtain that v is a distributional solution of

$ ' ' ' & ' ' ' % ´Bt vpt, xq ´σ2 2 ∆vpt, xq `H ``ş Ω pptq ˘Dvpt, xq şΩ pptq " f pxq ş Ω pptq , in Q T , vpT, xq " ´ 1 Ω pxq ş Ω ppT, yqdy , in Ω, (67)
i.e. that for all smooth function ψ :

Q T Ñ R with compact support in p0, T s ˆΩ, ż Q T vpt, xqB t ψpt, xqdxdt ` ş Ω ψpT, xqdx ş Ω ppT, xqdx `σ2 2 ż Q T Dvpt, xq ¨Dψpt, xqdxdt `żQ T H ``ş Ω pptq ˘Dvpt, xq şΩ pptq ψpt, xqdxdt " ż Q T f pxq ş Ω pptq ψpt, xqdxdt.
From this and (65), we deduce that

ż Q T upt, xqB t ψpt, xqdxdt ` ş Ω ψpT, xqdx ş Ω ppT, xqdx `σ2 2 ż Q T Dupt, xq ¨Dψpt, xqdxdt `żQ T H ``ş Ω pptq ˘Dupt, xq şΩ pptq ψpt, xqdxdt " ż Q T f pxq ş Ω pptq ψpt, xqdxdt ´żQ T ψpt, xqdγptqdx, (68)
which means that u P L 2 p0, T ; H 1 0 pΩqq is a distributional solution of (62), which furthermore satisfies (60). Finally ( 63) is a consequence of the weak convergence of p M pT q to ppT q in L 1 pΩq, the almost everywhere convergence of p M , the almost everywhere convergence of t Þ Ñ ş Ω p M ptq to t Þ Ñ ş Ω pptq, the almost everywhere convergence of p M |b M | 2 to p|b| 2 , using Fatou lemma when necessary.

6 Numerical method for the finite horizon problem 6.1 Numerical method Discretization.

We will use finite differences. To save notations, we focus on one-dimensional problems, but the same method can be generalized e.g. to dimension 2. Let us take Ω " r0, x max s Ă R for some x max ą 0. Let N T be a positive integer and let h P p0, x max q be such that x max {h is an integer. Let us denote ∆t " T {N T and N h " x max {h. We consider uniform grids on r0, T s and r0, x max s with respectively pN T `1q and pN h `1q points. Set t n " n∆t, and x i " i h for pn, iq P t0, . . . , N T u ˆt0, . . . , N h u. The values of u and p at px i , t n q are respectively approximated by U n i and P n i , for each pn, iq. We will use the notations U i " pU n i q n"0,...,N T , U n " pU n i q i"0,...,N h , and likewise for P . In the sequel, following [START_REF] Achdou | Finite difference methods for mean field games[END_REF] in the context of MFG, we consider a numerical Hamiltonian corresponding to q H defined by [START_REF] Davis | Algorithm 832: UMFPACK V4.3-an unsymmetric-pattern multifrontal method[END_REF]. More precisely, we consider r H : R ˆR˚ˆR ˆR Ñ R, px, µ, ξ 1 , ξ 2 q Þ Ñ r Hpx, µ, ξ 1 , ξ 2 q satisfying the following conditions:

• Monotonicity: r H is non-increasing with respect to ξ 1 and nondecreasing with respect to ξ 2 .

• Consistency: r Hpx, µ, ξ, ξq " q Hpx, µ, ξq.

• Differentiability: r H is of class C 1 in all variables.

• Convexity: for all x P R, µ P R ˚, pξ 1 , ξ 2 q Þ Ñ r Hpx, µ, ξ 1 , ξ 2 q is convex.

Remark 6.1. For example, if Lpx, bq " f pxq `1 2 |b| 2 , then q H is given in Remark 2.9 and for the numerical Hamiltonian, one can take

r Hpx, µ, ξ 1 , ξ 2 q " # 1 2 µ `pξ 1 q 2 `pξ 2 q 2 ˘´fpxq µ , if µ `pξ 1 q 2 `pξ 2 q 2 ˘1{2 ď M, 1 2 M `pξ 1 q 2 `pξ 2 q 2 ˘1{2 ´fpxq µ ´M2 2µ , otherwise, where x P R d , ξ 1 , ξ 2 P R, µ P R ˚.
We also consider discrete versions of F, G that we will be noted r F , r G : R N h `1 Ñ R N h `1. A typical set of assumptions that can be made on r ϕ P t r F , r Gu is that

• r ϕ is Lipschitz continuous from R N h `1 endowed with the discrete L 2 norm : }V } 2 " `h ř i V 2 i
˘1 2 to itself, with a Lipschitz constant independent of h • Let I h be the piecewise linear interpolation at the grid nodes. There exists a continuous and bounded function ω : R `Ñ R `such that ωp0q " 0 and for all p P L 2 pΩq, for all sequences pP phq q h , › › ›ϕrps ´Ih ´r ϕrP phq s

¯› › › L 2 pΩq
ď ω ´}p ´Ih P phq } L 2 pΩq ¯.

We introduce the finite difference operators: for W P R N T `1, pD t W q n " 1 ∆t pW n`1 ´W n q, n P t0, . . . N T ´1u, and for W P R N h `1, pD `W q i " 1 h pW i`1 ´Wi q, i P t0, . . . N h ´1u, p∆ h W q i " ´1 h 2 p2W i ´Wi`1 ´Wi´1 q , i P t1, . . . N h ´1u, r∇ h W s i " `pD `W q i , pD `W q i´1 ˘T , i P t0, . . . N h ´1u.

Discrete HJB equation.

To alleviate the notations, let us introduce, for P, U P R N h `1,

FpP, U q i " ´ÿ k hP k r H µ ˜xk , ÿ j hP j , DU k ¸`r F " P ř j hPj ı px i q ř j hP j ´ÿ k hP k r F " P ř j hPj ı px k q ´řj hP j ¯2
and

GpP q i " ´ ř j hP j `1 ř j hP j r G « P ř j hP j ff px i q ´řk P k r G " P ř j hPj ı px k qh ´řj hP j ¯2 .
Then, for the HJB equation ( 15), we consider the following finite difference scheme:

´pD t U i q n ´σ2 2 p∆ h U n q i `r H ˜xi , ÿ j hP n`1 j , r∇ h U n s i ¸" FpP n`1 , U n q i (69)
i P t1, . . . , N h ´1u, n P t0, . . . , N T ´1u,

U n i " 0, i P t0, N h u, n P t0, . . . , N T ´1u, (70)
U N T i " GpP N T q i , i P t0, . . . , N h u. (71)
Discrete FP equation. To define a discretization of the FP equation, we consider the weak form of (12 Equation (76) can be solved using backward time marching with Newton method at each time step. Equation (77) can be solved using forward time marching and solving a linear system of equations at each time step.

In our numerical implementation, in order to solve the linear systems associated to (76)-(77), we have used the open source library UMFPACK [START_REF] Davis | Algorithm 832: UMFPACK V4.3-an unsymmetric-pattern multifrontal method[END_REF] which contains an Unsymmetric MultiFrontal method for solving linear systems.

Numerical results

6.2.1 Case 1: a one dimensional problem with T " 0.2

We present here results in dimension 1 with a small time horizon: we take Ω " p0, 1q and T " 0.2. The cost functional is given by ( 6

The graphs of p 0 and g T are shown in Figure 1. The evolution of p and u is shown in Figure 2 for σ " 0.8. Since the goal is to minimize the cost functional, the mass moves towards the minimum of g T . The numerical convergence for this test case is illustrated by Figure 3. We have stopped the iterative procedure when the normalized 2 distance between two successive approximations of p and u were both smaller than 10 ´6, where the (time-space) normalized 2 norm is defined, for a vector V P R pN h `1qˆpN T `1q , by

}V } 2 " ˜h ∆t N T ÿ n"0 N h ÿ i"0 |V n i | 2 ¸1{2 .
For the results presented in Figures 123, we have used h " 5 ˆ10 ´4 and ∆t " 2 ˆ10 ´4. 6.2.2 Case 2: same problem with T " 2.0 Figure 4 displays the evolution of u and p with the same cost functions and initial distribution as above but with a much larger time horizon. The mass decays to 0 at an exponential rate but has a certain structure even at the final time, as shown in Figure 5. The larger T , the smaller ş Ω ppT, xqdx and hence the larger the final condition (78). This leads to numerical difficulties for very large time horizons. A possible approach is to rely on a continuation method, that is, to use the solution for a shorter time horizon to initialize the iterative procedure for a longer time horizon. However this method is time consuming since it requires solving a number of intermediate problems. We will come back to this matter in Section 7.3. For the results presented in Figures 45, we have used h " 5 ˆ10 ´4 and ∆t " 2 ˆ10 ´4.

Results in 2D

The numerical method described in Section 6 can be applied beyond dimension 1. Here, we provide results for two test cases in dimension 2. We take Ω " p0, 1q 2 and T " 0.2. For the results presented in this paragraph, we have used h " 1{80 " 0.0125 in both dimensions of space, and ∆t " 5 ˆ10 ´3. Case 3: First test case in dimension 2. The running cost is quadratic in the control. The initial distribution and the terminal cost are displayed in Figure 6: at the beginning the mass is concentrated around the center of the domain, but the final cost discourages to stay there. The evolution of the density is shown in Figure 7.

Case 4: Second test case in dimension 2. The running cost is quadratic in the control. The initial distribution and the terminal cost are displayed in Figure 8: at the beginning the mass is concentrated around the center of the domain, but the final cost attracts it towards two points. The evolution of the density is shown in Figure 9.

Numerical results

Stationary solution

We consider again the one-dimensional test case discussed in Section 6.2, except that we consider a running cost instead of a final cost. This example will be referred to as Case 5 in what follows. More precisely, the cost functional is given by (39), with Lpx, bq " nd " 0, for x P R, b P R, p P L 2 pΩ; Rq. The solution pp, uq is shown in Figure 10. It is compared with the solution at time T {2 of the problem with a time horizon T , rescaled by the remaining mass (namely, u and p are respectively multiplied and divided by ş Ω ppT {2q). The numerical convergence of the iterative procedure is illustrated by Figure 11. The stopping criterion is based on the normalized 2 distance between two successive approximations of p and u, where the normalized 2 norm is defined, for a vector V P R N h `1, by

}V } 2 " ˜h N h ÿ i"0 |V i | 2 ¸1{2 .
For the results presented in Figures 10-11, we have used h " 10 ´3.

Relation between non-stationary and stationary solutions

Figure 12 displays the evolution of (spatial) normalized 2 distance between the solution to the stationary problem and the solution to the time-dependent problem, with respect to time. The graphs are shown for several values of the time horizon T . One can see that, provided T is large enough, the distance tends to 0 for t bounded away from 0 and T . This is reminiscent of the so-called turnpike phenomenon: starting from the initial condition, the solution quickly evolves to a stationary regime, which it leaves around the end of the time interval in order to satisfy the final condition.

Alternative method for the finite time horizon problem

As described above (see Section 6.2.1), the fact that the total mass tends to 0 as time increases leads to numerical difficulties for large time horizons. However, the long time behavior suggests that a properly scaled version of the time dependent version should be more amenable to numerical treatment. As evidenced by

  bpt, X b t qq1 tτ b ątu ‰ Ppτ b ą tq " ş Ω Lpx, bpt, xqqp b pt, xqdx ş Ω p b ptq .Similarly, the final cost can be written asE " gpX b T q | τ b ą T ‰ " ş Ω gpxqp b pT, xqdx ş Ω p b pT q .Overall, the problem can be recast as the following deterministic control problem driven by a KFP equation:minimize ż T 0 ş Ω Lpx, bpt, xqqp b pt, xqdx ş Ω p b ptq dt `şΩ gpxqp b pT, xqdx ş Ω p b pT q subject to

Theorem 2 . 2 .

 22 For each b P B, the unique weak solution p b to (5) is continuous in Q T and positive in Q T . The cost Jpbq given by (6) is well defined for each b P B. Moreover, for every M ą 0, there exists b opt P B M which minimizes J over B M .

Theorem 3 . 5 (

 35 Theorem 2.4 in[START_REF] Berestycki | The principal eigenvalue and maximum principle for second-order elliptic operators in general domains[END_REF]). Consider a bounded domain Ω of R d , b P L 8 pΩ; R d q and a compact set K contained in Ω such that |K| ě δ 0 ą 0 (where |K| denotes the Lebesgue measure of K). Let λ b pΩq and λ b pΩzKq be the principal eigenvalues associated respectively with Ω and ΩzK. Then λ b pΩzKq ´λb pΩq ě δ ą 0,

Lemma 3 . 7 .

 37 Let the sequence b n tend to b in L 8 pΩ; R d q weak *. Let λ n be the principal eigenvalue of L n (related to b n ). Let p bn be the principal eigenvector of L n (non negative and normalized in L 1 norm). The following convergence results hold: λ n tends to λ b and p bn tends to p b in L 2 pΩq strongly and in H 1 0 pΩq weakly as n Ñ 8.

  ż Ω p bpxqLpx, bpxqqdx ď lim inf mÑ8 ż Ω p bm pxqLpx, bm pxqqdx and Φrp bs " lim mÑ8 Φrp bm s. Therefore, b is a minimizer of (38) subject to (40)-(41).

Ω

  pDu ¨Dδp b `δp bb ¨Du ´λb uδp bqdx `żΩ p bL bpx, bq ¨δb dx ` δλ b. Since ş Ω p be b ą 0, there exists a unique u b satisfying (46) and ż Ω p bu bdx " ´ .

3 .

 3 A positive Radon measure γ on r0, T s, such that, setting bpt, xq " ˆżΩ pptq ˙Dupt, xq 1.

Figure 1 :

 1 Figure 1: Case 1: data of the problem: initial density p 0 (left) and terminal cost g T (right).

Figure 2 : 2 Figure 3 :

 223 Figure 2: Case 1: evolution of p (left) and u (right), with T " 0.2.

Figure 4 :

 4 Figure 4: Case 2: evolution of p (left) and u (right), with T " 2.0.

Figure 5 :

 5 Figure 5: Case 2: evolution of t Þ Ñ ş Ω ppt, xqdx (left), and graph of ppT, ¨q at final time T " 2.0 (right).

Figure 6 :

 6 Figure 6: Case 3: data of the problem: initial density p 0 (left) and terminal cost g T (right).

Figure 7 :

 7 Figure 7: Case 3: evolution of the density.

Figure 8 :

 8 Figure 8: Case 4: Second test case in dimension 2. Data of the problem: initial density p 0 (left) and terminal cost g T (right).

Figure 9 :

 9 Figure 9: Case 4: Second test case in dimension 2. Evolution of the density.

Figure 10 :

 10 Figure 10: Case 5: Solution to the stationary problem: p (left), u (right) and λ « 3.15.

Figure 11 :

 11 Figure 11: Case 5: Numerical convergence: normalized 2 distance between two iterations, for p (top left) and u (top right), and evolution of ş Ω pu (bottom) with respect to the number of iterations.

Figure 12 :

 12 Figure 12: Case 5: Comparison between the stationary and non-stationary solutions for 3 different time horizons (T " 0.5, 1, 2): normalized 2 distance between stationary and non-stationary p (left) and u (right).

  In this section M is fixed. Later on, when necessary, we will change the notation to make the dependency of the Hamiltonian on M explicit. Theorem 2.6. Let b opt P B M be a minimizer of J over B M . Then for almost every pt, xq P Q T , b opt pt, xq achieves the maximum in[START_REF] Boccardo | Existence de solutions non bornées pour certaines équations quasi-linéaires[END_REF] with ξ replaced by `şΩ ppt, yqdy ˘Dupt, xq, namely,

	Remark 2.5. b opt pt, xq " H ξ ˆx,	ˆżΩ	pptq ˙Dupt, xq	˙,
	where pp, uq solves the forward-backward PDE system
	$ ' ' ' ' ' ' ' &	´Bt upt, xq	´σ2 2	∆upt, xq	`H `x, `şΩ pptq ˘Dupt, xq pptq	şΩ " " ppt,¨q F ş Ω pptq
	'
	'
	'
	'
	'
	'
	'
	%
							in the form of a forward-
	backward PDE system. Let us introduce the Hamiltonian H : Ω ˆRd Ñ R,
					Hpx, ξq "	max bPR d ; |b|ďM	pξ ¨b ´Lpx, bqq .	(10)

  Solving the finite difference equations.

						).
	It involves, for a smooth w,
		´żΩ	div ˆppt, ¨q q H ξ ˆ¨,	ż	pptq, ∇upt,	¨q˙˙p xqwpt, xqdx
						Ω
	"	ż	ppt, xq q H ξ ˆx,	ż
		Ω

Ω pptq, ∇upt, xq ˙¨∇wpt, xqdx,

Acknowledgment. The research of the first author was partially supported by the ANR (Agence Nationale de la Recherche) through MFG project ANR-16-CE40-0015-01.

(52)

Note that by contrast with [START_REF] Boccardo | Existence de solutions non bornées pour certaines équations quasi-linéaires[END_REF], the Hamiltonians in (52) do not involve f because it is more convenient for what follows.

The optimality conditions for (50) consist of a forward-backward system of PDEs:

and γ M : r0, T s Ñ R is defined by

Remark 5.1. For any smooth function φ : r0, T s Ñ R compactly supported in p0, T q, we test (53) by p M φ and (54) by u M φ, subtract the resulting identities and integrate by part: we obtain that

We wish to study the behaviour of a family of solutions pb M , p M , u M q M as M Ñ `8.

Theorem 5.2. There exist 1. u P L 2 p0, T ; H 1 0 pΩqq X L 8 pQ T q 2. p P C 0 pr0, T s; L 1 pΩqq X L 8 p0, T ; LlogLpΩqq X L r pQ T q X W 1,q p0, T ; W ´1,q pΩqq for all

where we used integration by parts and the Dirichlet boundary condition.

This leads us to introduce the following discrete operator (see [START_REF] Achdou | Finite difference methods for mean field games[END_REF] for more details), for µ P R, U, P P R N h `1, and for i P t1, . . . , N h ´1u, B pµ,U q i pP q " 1 h

Then, for the discrete version of the FP equation ( 12) we consider the following finite difference equation:

i P t1, . . . , N h ´1u, n P t0, . . . , N T ´1u,

Numerical method.

An interesting option for solving ( 11)-( 12) is to use an Augmented Lagrangian method as was done in [START_REF] Achdou | Mean field type control with congestion II : An augmented Lagrangian method[END_REF] for mean field type control problems, but we have not done this yet. One could also try to solve directly the whole system ( 11)-( 12) using for example Newton method, see [START_REF] Achdou | Finite difference methods for mean field games[END_REF] in the context of MFG. However, equations ( 69) and ( 73) involve some non-local terms in P and U , and hence one would need to invert full matrices, which would probably be very time consuming. For this reason we instead use the following iterative method, where the solution computed at a given iteration is used to compute the non-local terms involved in the next iteration:

1. Start with a guess pP p0q , U p0q q; set k Ð 0.

2. Repeat the following steps (a) Compute U pk`1q , solution to the following modified version of (69)-(71),

, r∇ h U n s i ¸" FpP pkq,n`1 , U pkq,n q i , i P t1, . . . , N h ´1u, n P t0, . . . , N T ´1u,

(b) Compute P pk`1q , solution to the following modified version of (73)-(75),

,n q i pP n`1 q " 0, i P t1, . . . , N h ´1u, n P t0, . . . , N T ´1u,

(c) If ||P pk`1q ´P pkq || 2 and ||U pk`1q ´U pkq || 2 are small enough, stop. Otherwise set k Ð k `1 and continue.

To ensure convergence, it is sometimes helpful to use relaxation in the update rule, that is, denoting by Ũ pk`1q the solution of (76), we set U pk`1q " p1 ´θqU pkq `θ Ũ pk`1q , for some θ P p0, 1q, and similarly for P pk`1q .

7 Numerical method for the stationary problem

Numerical method

To solve ( 16)-( 17), we first notice that we can replace this system by the system composed of ( 16) which we rewrite for convenience

2 ∆q p ´div pq pH ξ p¨, Dq uqq " λq p, in Ω, q p " 0 on BΩ,

(79) and the following modified HJB equation:

2 ∆q u `Hp¨, Dq uq " λq u `F rq ps `c1 1 Ω ´ ´żΩ q uq p, in Ω,

´`ş Ω q p pL p¨, H ξ p¨, Dq uqq `F rq psq dx ˘.

(80) Indeed ( 16)-( 17) is equivalent to (79)-( 80): it is straightforward to see that ( 16)-( 17) implies ( 79)-( 80); conversely multiplying the HJB equation ( 80) by q p and the Fokker-Planck equation ( 79) by q u, integrating on Ω, summing the resulting equations and using the definition of c 1 , we see that most of the terms cancel and we obtain ż Ω q upxqq ppxqdx " ´ and finally ( 16)-( 17). This allows us to use the same kind of iterative algorithm as for the non stationary problem (see Section 6), except that, in order to solve (79), the linear solver for the FP equation is replaced by an inverse power method. More precisely, we use the following iterative method, where the solution computed at a given iteration is used to compute the non-local terms involved in the next iteration. We will use the notations r H and F introduced in Section 6.1 for the finite time horizon problem. The iterative procedure we consider is:

1. Start with a guess p q P p0q , q U p0q , λ p0q q; set k Ð 0.

2. Repeat the following steps (a) Compute q U pk`1q , solution to the following modified version of (80),

`Fp q P pkq , q U pkq q i ´ ´h ř j q U pkq j q P pkq j , i P t1, . . . , N h ´1u,

(b) Compute p q P pk`1q , λ pk`1q q, solution to the following modified version of (79),

2 p∆ h q P q i ´Bp1, q U pk`1q q i p q P q " λ q P i , i P t1, . . . , N h ´1u, q P i " 0, i P t0, N h u, q P i ě 0, i P t1, . . . , N h ´1u, ř j h q P j " 1.

(82) (c) If || q P pk`1q ´q P pkq || 2 and || q U pk`1q ´q U pkq || 2 are small enough, stop. Otherwise set k Ð k `1 and continue.

To solve (81), here again we use UMFPACK. For the second step, namely finding the eigenvalue and eigenvector solving (82), our implementation relies on the ARPACK routines [START_REF] Lehoucq | ARPACK users' guide, volume 6 of Software, Environments, and Tools[END_REF], which are based on the Implicitly Restarted Arnoldi Method.

Figure 12, the proper scaling should come from the total remaining mass, but this is unknown until we solve the non-stationary problem. Let us instead rescale the functions using the ergodic constant coming from the non-stationary problem. More precisely, for γ ą 0, let pp pT,γq , u pT,γq q solve the following scaled system

and c 1 , c 2 : r0, T s Ñ R are defined by [START_REF] Carmona | Probabilistic theory of mean field games with applications[END_REF].

In order to solve the original system ( 11)-( 12), we propose the following method:

1. Solve the non-stationary problem ( 16)-( 17) and obtain pp, ũ, λq.

2. Solve the scaled system with parameter γ " λ and obtain pp pT,λq , u pT,λq q.

3. Un-scale the solution by letting u pT q pt, xq " e λt u pT,λq pt, xq and p pT q pt, xq " e ´λt p pT,λq pt, xq.

One can check that pp pT q , u pT q q indeed solves ( 11)-( 12).