Maximal $L^q$-regularity for parabolic Hamilton-Jacobi equations and applications to Mean Field Games

In this paper we investigate maximal $L^q$-regularity for time-dependent viscous Hamilton-Jacobi equations with unbounded right-hand side and superlinear growth in the gradient. Our approach is based on the interplay between new integral and H\"older estimates, interpolation inequalities, and parabolic regularity for linear equations. These estimates are obtained via a duality method \`a la Evans. This sheds new light on a parabolic counterpart of a conjecture by P.-L. Lions on maximal regularity for Hamilton-Jacobi equations, recently addressed in the stationary framework by the authors. Finally, applications to the existence problem of classical solutions to Mean Field Games systems with unbounded local couplings are provided.


Introduction
The purpose of this paper is two-fold. First, to establish maximal L q -regularity for parabolic Hamilton-Jacobi equations with unbounded right-hand side of the form ∂ t u(x, t) − ∆u(x, t) + H(x, Du(x, t)) = f (x, t) in Q T = T d × (0, T) , where H has superlinear growth in the gradient variable and f ∈ L q (Q T ) for some q > 1. Second, to apply these regularity results to prove the existence of classical solutions for a large class of second order Mean Field Games systems with local coupling, i.e.
Maximal regularity for (HJ). The problem of maximal L q -regularity for (HJ) amounts to show that bounds on the right-hand side f in L q (Q T ) imply bounds on individual terms ∂ t u, ∆u and H(x, Du) in L q (Q T ) (see Theorem 1.1 below for a more precise statement). This problem has been proposed for stationary Hamilton-Jacobi equations by P.-L. Lions in a series of seminars (see e.g. [42,41]). Under the assumption that H(x, p) = |p| γ , γ > 1, he conjectured that maximal regularity holds provided that q is above the threshold d(γ − 1)/γ. The conjecture has been proved recently in [21] via a refined Bernstein method, which unfortunately breaks down in the parabolic setting. Here, we are able to obtain parabolic maximal regularity via different methods, assuming that q is above a certain (parabolic) threshold, that is in the regime of sub-natural growth, that is for γ < 2. For γ ≥ 2, we obtain the larger threshold for every x ∈ T d , p ∈ R d . Concerning the case γ ≥ 2, we will further assume some additional regularity in the x-variable, i.e for α ∈ (0, 1) to be specified, for all x, ξ ∈ T d and p ∈ R d . A typical example of H satisfying (H) is If h ∈ C α (T d ), this Hamiltonian will satisfy also (H α ).
Hoping to help the reader to have a clearer picture, we sketch known and new regularity regimes as γ and q vary in Figure 1. Theorem 1.1. Assume that (H) holds, and (H α ) also when γ ≥ 2 (with α as in Theorem 1.2 below). Let u ∈ W 2,1 q (Q T ) be a strong solution to (HJ) and assume that for some K > 0 then, there exists a constant C > 0 depending on K, q, d, C H , T such that The strategy of the proof is based on the following procedure. By maximal regularity for linear equations, one has D 2 u L q H(Du) L q + f L q + u 0 W 2−2/q, q Du γ L γ q + f L q + u 0 W 2−2/q, q . Then, one looks for a suitable norm |||·||| so that a Gagliardo-Nirenberg type interpolation inequality of the form Du γ L γ q D 2 u γθ L q |||u||| γ(1−θ) with γθ < 1 holds. Combining the two inequalities, maximal regularity is achieved whenever it is possible to produce bounds on |||u|||. This way of producing estimates for nonlinear problems, using linear estimates and interpolation inequalities, goes back to the works of Amann and Crandall [2]. When γ = 2, a good choice for |||·||| is the L ∞ norm, which is easily controlled by f ∞ (ABP type estimates allow to reach f d+1 , as in [58]). Here, we start with the observation that when γ < 2, the optimal choice is a suitable L p norm, while for γ > 2, a Hölder bound on u is needed. A crucial step in this work is the derivation of such estimates. These are obtained by duality arguments, inspired by [28,20]. The main idea behind duality is to shift the attention from (HJ) to the formal adjoint of its linearization, which is a Fokker-Planck equation. Crossed information on u and the solution ρ to the "dual" equation allow to retrieve estimates on ρ, that are then transferred back to u. A more detailed heuristic explanation of this method can be found in [20] (see also references therein). Note that [20] is devoted to Lipschitz regularity of u, while here we investigate Hölder regularity, which requires a further inspection of the regularity of ρ at the level of Nikol'skii spaces.
We underline that our results on Hölder regularity in the super-quadratic case are new in the following sense. Recall that Hölder bounds have been obtained in [13,60] when q > 1 + d/γ, while here we assume the stronger requirement q > (d + 2)/γ (which is the natural one for maximal regularity). In [13,60] sign assumptions on f are in force, and explicit Hölder exponents are not provided. Here, we do not require any assumption on the sign of f , and produce explicit Hölder exponents. The statement is as follows.
We mention that the proof of the Hölder bounds could be localized in time, thus assuming merely u 0 ∈ C(T d ). The constant C will then depend on u ∞ , see Remark 3.7.
In the sub-quadratic regime γ < 2, the existence and uniqueness of weak solutions can be obtained up to the critical integrability exponent q = (d + 2) γ−1 γ [44]. From the viewpoint of maximal regularity, this endpoint situation is a bit more delicate than the one treated in Theorem 1.1, which concerns q > (d+2) γ−1 γ . Indeed, the heuristic procedure previously discussed yields which is meaningful only if u L p is small. We circumvent this issue by shifting the analysis from u to u −u k , where u k is the solution to a suitable regularized problem. Crucial stability estimates on u −u k L p then lead to the second main maximal regularity result of the paper.

Theorem 1.3.
Assume that (H) holds, and 1 + 2 d+2 < γ < 2. Let u ∈ W 2,1 q (Q T ) be a strong solution to (HJ), and We stress that in this limiting case the dependance of the constant C with respect to f is not just through its norm, as in the previous Theorem 1.1, but on subtler properties; see Remark 4.2 below for additional details.
Further comments concerning the thresholds for q in Theorems 1.1 and 1.3 are now in order. When q < (d + 2)/γ we believe that maximal regularity results are false in general. This would be in line with known results for the associated stationary problem [21], for which maximal regularity does not hold when q < d/γ . We also believe that the endpoint q = d+2 2 for quadratic problems γ = 2 could be treated by our methods, using more refined interpolation estimates in Orlicz spaces (or the analysis of a linear equation obtained via the Hopf-Cole transformation). Regarding the super-quadratic regime γ > 2, we do not know whether our assumption q > (d + 2)(γ − 1)/2 is optimal or not.
It is finally worth noticing that our results, and in particular the Hölder estimates in Theorem 1.2, apply to equations with repulsive gradient term (e.g. Kadar-Parisi-Zhang type PDEs), i.e.
with G satisfying (H). In other words, the sign in front of H and f does not matter in (HJ), since it is sufficient to observe that u(x, t) = −v(x, t) solves (HJ) with H(x, p) = G(x, −p), which satisfies (H) too. We refer to [4,5] and the references therein for further discussions on these equations.

Mean Field Games.
Armed with maximal regularity results for Hamilton-Jacobi equations, we describe our contributions on Mean Field Games (MFG) systems of the form (MFG). These arise in the MFG framework, which is a class of methods inspired by ideas in statistical physics to study differential games with a population of infinitely many indistinguishable players [39,38]. The PDE system (MFG) appears naturally when the running cost of a single player depends on the population's density in a pointwise manner. This results in the nonlinear coupling term g(m(x, t)) between the Hamilton-Jacobi and the Fokker-Planck equation in (MFG). When g(·) is an unbounded function, regularity of solutions is still a challenging problem that has not been solved in full generality.
In the so-called monotone case, that is when g(·) is increasing, there are basically no restrictions on the growth of f if one looks for solutions in some weak sense [11,53], or for classical solutions when the time horizon T is small [3,18]. Classical solutions for arbitrary T can be obtained by requiring a mild behaviour of g(m) as m → ∞, or a mild behaviour of H(p), i.e. γ ≤ 1 + 1/(d + 1) (as suggested in [41]). We deal here with the first situation. For the model problem g(m) = m r , an upper bound on r depending on γ and the space dimension d is usually required. Concerning the subquadratic case γ < 2, a main reference is [30] (see also [6]), while [32] explores the superquadratic case 2 < γ < 3.
We will assume here that f : T d × [0, ∞) → R is of class C 1 , and that there exist r > 0 and C g > 0 such that This implies that g(·) is monotone increasing, and bounded from above and below by power-like functions of type m r . As for the Hamiltonian, we will further assume some additional smoothness and uniform convexity in the p-variable, having always in mind a power-like growth |p| γ : for every p ∈ R d , M ∈ Sym d Our main result on this class of problems reads as follows.
It is understood here that there are no restrictions on r when d = 1, 2. To our knowledge, the results of this manuscript extend known classical regularity regimes in the sub-quadratic case. Note that Regarding the super-quadratic case, we actually obtain the same restriction r < 2 d(γ−1)−2 as in [32], but we are able to cover the full interval γ ∈ [2, ∞), thus unlocking some smoothness regimes for γ > 3. Still, one may conjecture that, for all d ≥ 1, no restrictions on r should be required to get classical solutions (as for the purely quadratic case H(p) = |p| 2 investigated in [12, Theorem 4.1]), but this remains an open question.
In the non-monotone framework, the assumption of increasing monotonicity of g(·) is dropped. Conversely, one may even consider a g(·) which is strictly decreasing (focusing case), and model aggregation phenomena as in [14,16,22]. In this framework, a control on the growth of g is structurally needed. For stationary problems, it was shown in [15] that when g(m) = −m r , no solutions on R d exist when r > γ d−γ , and even for r ≥ γ d some further assumptions might be needed for solutions to be obtained. On the other hand, when r < γ d , existence of weak solutions is shown in [22], but their classical regularity is proven under much stronger hypotheses, that impose γ < 2. Now, we suppose With respect to the previous assumptions (M + ) in the monotone case, g need not be monotone increasing; in contrast, it can be strictly decreasing, for example g(m) = −m r . We have the following In other words, we prove that solutions in the focusing case obtained in [22] are always classical in the subquadratic regime. In the superquadratic regime, a new class of smooth solutions is obtained. We stress that the threshold γ d is crucial, as one may have non-existence of solutions for large time-horizon T when r ≥ γ d [17]. For a summary sketch of the existence regimes that we obtain here, see Figure 2. Theorems 1.4 and 1.5 are proven as follows. First, we use structural estimates to deduce some a priori bounds on g(m) L q . These are mainly second-order type estimates in the monotone framework, and first-order estimates in the non-monotone one (plus some additional interpolation procedure). The assumptions on the growth r of g(·) guarantee that q is large enough to apply our maximal regularity results for the HJ equation. Then, once the crucial a priori bounds on u in W 2,1 q are established, a bootstrap procedure allows to get estimates up to second order derivatives in Hölder spaces, and existence follows via standard methods.
Outline. In Sections 2.1 and 2.2 we introduce the functional setting that is used throughout the paper, while Section 2.3 is devoted to Sobolev regularity results for Fokker-Planck equations. Section 3 comprehends the proof of a priori integral and Hölder estimates for Hamilton-Jacobi equations with unbounded ingredients. In Section 4, Theorem 1.1 on maximal regularity for (HJ) is proven, while in Section 5 the existence Theorems 1.4 and 1.5 for (MFG) are proven. work has been partially supported by the Fondazione CaRiPaRo Project "Nonlinear Partial Differential Equations: Asymptotic Problems and Mean-Field Games". The second-named author wishes to thank Prof. Alessandra Lunardi for discussions about Lemma A.1.

Fractional spaces of periodic functions
Let µ ∈ (0, 1) and and it is endowed with the equivalent norm is known as Nikol'skii space [49] and the above norm is interpreted as usual as Then, we define the Nikol'skii class We mention that in view of [63, Corollary 2 p. 143], Nikol'skii spaces N µ, p (T d ) can be endowed with equivalent norms of the form where k, r ∈ Z are such that 0 ≤ k < µ and r > s − k, where ∆ r h are the increments of order r and step h, see [40]. Moreover, these spaces can be characterized via real interpolation: for m ∈ N, p, q ∈ [1, ∞] and θ ∈ (0, 1) we have1 , with equivalence of the respective norms, see e.g. [43], [40,Theorem 17.24].

Space-time anisotropic spaces
. We will also use the notation Q t 2 : x u ∈ L p (Q) for all multi-indices β and r such that | β| + 2r ≤ 2, endowed with the norm The space W 1,0 p (Q) is defined similarly, and is endowed with the norm .
We define the space H 1 p (Q) as the space of functions u ∈ W 1,0 ; X) and L q (t 1 , t 2 ; X) the usual spaces of continuous, Hölder and Lebesgue functions respectively with values in a Banach space X, we have the following isomorphisms: and the latter is known to be continuously embedded into C([t 1 , t 2 ]; L 2 (T d )) (see, e.g., [25, Theorem XVIII.2.1]). Sometimes, we will use the compact notation C(X) and L q (X).

Sobolev regularity of solutions to Fokker-Planck equations
In this section we collect a few existence and regularity properties of solutions to 1We denote here with (X 0 , X 1 ) θ, q the standard real interpolation space between Banach spaces X 0 , X 1 .
, and ρ τ be as in (5). Then, there exists a unique weak for all s > 0 and ϕ ∈ H 1 2 (Q s,τ ), and ρ(τ) = ρ τ in the L 2 -sense. Moreover, ρ is a.e. nonnegative on Q τ . Our analysis will be based on regularity properties of ρ in Sobolev spaces that depend on the integrability of |b| against ρ itself. Similar results already appeared in [22,9,47,53].

Proposition 2.2.
Let ρ be the nonnegative weak solution to (4) and 1 < σ < d + 2. Then, there exists where As we will se below, if σ > d+2 d+1 we obtain also Proof. The case σ < d+2 d+1 is covered by [20,Proposition 2.5], which is based on duality combined with maximal regularity arguments (following [47]). We focus here on the case σ > d+2 d+1 , and prove the theorem via a (standard) method from weak solutions, that does not exploit parabolic Caldèron-Zygmund regularity. This is possible because d+2 To simplify, we use ϕ := ρ β as a test function in the weak formulation of (4) integrating on Q t,τ := T d × (t, τ), while the argument can be made rigorous by testing against ϕ := (ρ + ) β , and then letting → 0. We have In particular, noting that Dρ We first pass to the supremum over t ∈ (0, τ) and, by means of [27, Proposition I. 3.1] we deduce We then have We then note that σ < d + 2 implies m > 2 and We apply Young's inequality to the first term on the right-hand side to conclude . Hence We thus conclude the estimate Finally, recalling that σ = d+2 d+3−m , we get the Sobolev estimate applying Hölder's inequality, using (9) and finally exploiting Young's inequality as .

Corollary 2.3. Let ρ be the nonnegative weak solution to
Proof. The first estimate follows by Proposition 2.2, applied with m = d+2 q (see also (8)), and the continuous embedding of The second estimate follows analogously, using also the embedding of

Bounds of solutions to HJ equations in Lebesgue and Hölder spaces
We derive in this section some preliminary bounds on solutions to Hamilton-Jacobi equations via duality methods. Since we are in the setting of maximal regularity and f ∈ L q (Q T ), we assume that u ∈ W 2,1 q (Q T )∩W 1,0 γq is a strong solution to (HJ), i.e. it solves the Hamilton-Jacobi equation almost everywhere. The initial datum u 0 is achieved in W 2−2/q,q (T d ), as suggested by the Lions-Peetre trace method in interpolation theory (cf Lemma A.1 below). Note that by classical embedding properties for Sobolev-Slobodeckij spaces, we have the following inclusions Since we will work under the assumptions in the statements of our results), and u ∈ L 2 (0, T; W 1,2 (T d )). Furthermore, Du(x, t)). Note finally that when q > (d + 2)/2, any solution is automatically continuous and u solves (HJ) in the weak sense used in [20]. This always happens in the superquadratic regime γ > 2. On the other hand, in the subquadratic case γ < 2, u is not necessarily continuous when (d + 2)/γ ≤ q < (d + 2)/2.
Remark 3.1. The assumption γ > 1 + 2 d + 2 guarantees that u ∈ L 2 (0, T; W 1,2 (T d )), so u has finite energy and it can be safely used as a test function for the dual equation (4). One can drop this requirement in all of the following statements, relaxing to so that q > d+2 γ > 1, and assuming a priori that u ∈ L 2 (0, T; W 1,2 (T d )) (which is always true for example when u is a classical solution, as in Section 5 on MFG). One could also invoke methods from renormalized solutions, to deal also with the case γ ≤ 1 + 1 d+1 but this is beyond the scopes of this paper.

Bounds in Lebesgue spaces
We obtain in this section bounds on the L p -norm of u, elaborating in particular the case d+2 γ ≤ q < d+2 2 , that is when γ < 2.
T k (s) = max − s, min{s, k} below will denote the truncation operator at level k > 0, u + = max{u, 0} and u − = (−u) + the positive and negative part of u respectively.
We start with some bounds on u + , that are obtained with no restrictions on q ≥ 1.

Lemma 3.2.
Assume that H is non-negative. Let u be a strong solution to (HJ) in W 2,1 q (Q T ), q ≥ 1. There exists a positive constant C 0 (depending on q, d, T) such that Note that W 2,1 q (Q T ) is embedded into C([0, T]; L p (T d )) (p as in the previous statement) so (11) has the form of an a priori bound.
Proof. We detail the proof in the case q < d+2 2 only. The case q > d+2 2 can be treated in an analogous way (essentially following [20,Section 3]).
For k > 0, let µ = µ k be the weak non-negative solution of the following backward problem Note that µ τ L p (T d ) ≤ 1. Since µ solves an equation of the form (4) with b ≡ 0, by Corollary 2.3 we have where C does not depend on k. Since u + is a weak subsolution to testing against µ, and testing the equation for µ against u we get that H(x, Du)µ dxdt.
We apply Hölder's inequality to the second term of the right-hand side of the above inequality, the assumption H ≥ 0 on the Hamiltonian, and the fact that the backward heat equation preserves the L p norm, i.e. µ(t) L p (T d ) ≤ 1 for all t ∈ [0, τ], to get, after sending k → ∞, the desired inequality.
We now proceed with some more delicate bounds on u − , under the restriction q > d+2 γ . We will use the following property of H: under the standing assumptions on H, the Lagrangian L(x, v) = sup p {ν · p − H(x, p)} satisfies for some C L > 0 (depending on C H ) Lemma 3.3. Assume that (H) holds. Let u be a strong solution to (HJ) in W 2,1 q (Q T ), q > d+2 γ . There exists a positive constant C (depending on C H , q, d, T) such that Proof. As before, we detail the case q < d+2 2 only. For k > 0, let ρ = ρ k be the weak non-negative solution of where C does not depend on k. Since u − is a weak subsolution to testing against ρ, and testing the equation for ρ against u − we get that On one hand, in view of (12) note that and on the other hand by Hölder's inequality Then, plugging (14) into (15) we obtain Since q > d+2 γ one can use Young's inequality, and the fact that ∫ ρ(t)dx ≤ ρ(t) L p (T d ) ≤ 1 to get the desired inequality (after passing to the limit k → ∞).

Corollary 3.4.
Assume that (H) holds. Let u be a strong solution to (HJ) in W 2,1 q (Q T ), q > d+2 γ , and assume that Then, there exists a positive constant C (depending on K, C H , q, d, T) such that Note that the previous result does not cover the critical case q = d+2 γ . The rest of the section is devoted to this endpoint situation, and more precise information on the stability of solutions in L p will be obtained. This will be crucial in the subsequent analysis of maximal regularity. It is worth noting that constants appearing in estimates below will not depend just on the norms f L q , u 0 L p , but on finer properties of f in L q and u 0 in L p ; see Remark 3.6.
Let Γ(x, t) be the fundamental solution of the heat equation on T d × (0, ∞). Consider, for k > 0, the solution u k to the problem with truncated / regularized data, i.e.
In other words, u k solves an Hamilton-Jacobi equation with L ∞ right-hand side and initial datum in C ∞ . Such an initial datum is actually u k (x, 0) = z(x, 1/k), where z is the solution to the heat equation and converges to u 0 in W 2−2/q,q (T d ) as k → ∞. The existence and uniqueness of a strong solution u k ∈ W 2,1 p (Q T ) for all p ≥ 1 can be obtained using for example results in [20]. We prove now estimates on u − u k .
Proof. Let w = u − u k . Note that w depends of course on k, but we will drop the subscript for simplicity. As before, we argue by duality, and estimate w + first. Fix τ ∈ (0, T] and ρ = ρ k be the weak nonnegative solution of Note that as in previous lemmas one should further truncate ρ(τ) to ensure the existence of ρ in an energy space (cf Proposition 2.1), and then pass to the limit, but we will omit this step for brevity.
Step 1: bounds on ∬ |D p H(x, Du k )| γ ρ. Testing (18) against ρ, and testing the equation for ρ against u k we get that Recall that for all t, Then, for h > 0 that will be chosen below, Similarly, for any h 0 > 0, by Young's inequality for convolutions and Hölder inequality, Finally, applying Lemma 3.2 to u k and noting that (u 0 Γ(1/k)) + ≤ u + 0 Γ(1/k) by the comparison principle, Pluggin the previous inequalities back in (20), and using the bounds from below on L, we then get for some C 1 depending only on T, d, γ. Finally, h and h 0 are chosen large enough so that Step 2: bounds on w + . Taking the difference between (18) and (HJ), by convexity of H(x, ·), w + is a weak subsolution of Testing against ρ, and testing the equation for ρ against w + gives Hence, using Corollary 2.3 and the estimate (21) we obtain which is "half" of the desired estimate.
To get an analogous estimate for w − L p (T d ) , which allows to conclude since w = u − u k , one can proceed as in Step 1 and 2, noting that w − satisfies To argue as before, it is sufficient to exploit properties of the dual problem Remark 3.6. Note that (19) directly yields the estimate thus extending (17) up to the critical case q = (d + 2)/γ . Let us focus on the the way C depends on f and u 0 . When q > (d + 2)/γ , C depends on f L q and u 0 L p only. At the endpoint q = (d + 2)/γ , C above (and similarly the constant in (19)) is proportional to C defined in (21), that depends in turn on and C 1 = C 1 (T, q, d) is as in Corollary 2.3. Thus, these constants remain bounded when f and u 0 vary in bounded and equi-integrable sets in L q (Q T ) and L p (T d ) respectively. This is completely in line with results in [44].

Bounds in Hölder spaces
We now proceed with bounds on u in Hölder spaces, which will be obtained in particular when γ ≥ 2, that is when q is necessarily greater than d+2 2 (and therefore u is continuous). Proof of Theorem 1.2. Step 1. Since we have the representation H(x, Du(x, t) (22) for any measurable Ξ : Q τ → R d such that L(·, Ξ(·, ·)) ∈ L r (Q τ ) and Ξ · Du ∈ L r (Q τ ), r > 1, and test function ϕ ∈ H 1 2 (Q τ )∩L r (Q τ ). The previous inequality becomes an equality if Ξ( Fix now any τ ∈ [0, T], and letx,ȳ ∈ T d be such that Let ρ τ be any smooth non-negative function satisfying Use now (22) with Ξ(x, t) = D p H(x, Du(x, t)) and ϕ = ρ, and u ∈ H 1 2 (Q T ) as a test function for the equation satisfied by ρ to get As before, plugging Ξ(x, t) = D p H(x − ξ, Du(x − ξ, t)) and ϕ =ρ into (22), and using u ∈ H 1 2 (Q T ) as a test function for the equation satisfied byρ yields Du(x, t)))ρ dxdt + ∬ Q τ fρ dxdt. (24) Taking the difference between (24) and (23) we obtain Step 2. To estimate the terms appearing in the right hand side of (25), we first derive bounds on ρ. We stress that constants C, C 1 , . . . below are not going to depend on τ and ρ τ . Rearranging (23) we have and by (12) and bounds on u ∞ of Proposition 3.4 we get Using Corollary 2.3, This provides a control on ∬ Q τ |D p H(Du)| γ ρ, and by means of Proposition 2.2, Step 3. First, recalling that As for the second term in (25), Du(x, t)).
Next, using (H α ), Finally, we apply the embeddings of Propositions A.3 and A.2 (with δ = q , p = d+2 d+3−γ and hence α = γ − d+2 q ) to get Plugging now all the estimates in (25) and using (26) we obtain It is now sufficient to recall that ρ τ can be any smooth non-negative function satisfying and we have the assertion.
Remark 3.7. Adding an additional time localization term in the previous procedure, as in [20], it is possible to obtain Hölder bounds that are independent of the initial datum, but just depend on the sup-norm of the solution. This indicates that the equation regularizes at Hölder scales, and weak solutions (in an appropriate sense) become instantaneously Hölder continuous at positive times.

Maximal L q -regularity
We start with a straightforward consequence of parabolic regularity results for linear equations.
for some positive constant C depending on q, d, C H .
Proof. The proof is an easy consequence of well-known Caldèron-Zygmund type maximal regularity results for heat equations with potential which satisfies the estimate (see [37,36], or [33] and the references therein) ).
To get (27), it is now sufficient to choose V = −H(x, Du) + f , and use the assumption (H).
We now proceed with our main result on maximal regularity for (HJ) when q > (d + 2)/γ .
We now consider the maximal regularity problem in the limiting case q = (d + 2)/γ . The scheme of the proof is similar to the one of Theorem 1.1, but requires an additional step involving solutions u k to the regularized problem (18), that is; for k > 0, Proof of Theorem 1.3. Let w = u − u k , k to be chosen. From Proposition 3.5, we have the existence of C depending on f , u 0 , C H , q, d, T such that The Gagliardo-Nirenberg inequality reads where C 2 depends on d, p, q. Thus, Note now that w solves a.e. on Q T and that, by assumptions on H, |D p H(x, p)| ≤ C H (|p| γ−1 + 1), so by Young's inequality where C 3 , C 4 depend on C H only. Then, by maximal regularity applied to the linear equation (34), where C 6 depends on C H , d, q. Plugging now this inequality into (33) yields Hence, we choosek large enough so that to get Since uk(0) is smooth and Tk( f ) ∈ L ∞ (Q T ), we can apply Theorem 1.1 to uk solving (18) to estimate Duk. Indeed, pick anyq > q. Then, in view of standard decay estimates for the heat equation (C 5 depends on d, q,q only, see e.g. [62,Chapter 15]). Therefore, by Theorem 1.1, Duk where Ck depends onk, , q, d, C H , T. Actually, Du can be proven to be bounded in L ∞ (Q T ), see [20]. It is now straightforward to conclude. Indeed, and the assertion follows by (36) and (37).
We claim that C appearing in the statement of Theorem 1.3 remains bounded when • f varies in a bounded and equi-integrable set F ⊂ L q (Q T ), and Indeed, in addition to , q, d, T, C H , the constant C crucially depends onk appearing in (35).
This is chosen in the proof large enough so that In turn, c = (2C 6 C γ 2 C γ−1 ) −1 is independent of f ∈ F and u 0 ∈ U 0 , since they vary in bounded and equiintegrable sets in L q (Q T ) and L p (T d ) respectively, cf. Remark 3.6. Note that by Sobolev embeddings, the closure of U 0 in L p (with p, q as above) is compact in L p (T d ), and hence weakly compact and L p -equi-integrable by the Dunford-Pettis theorem, see [10,Theorem 4.30].
Hence, we just need to verify that for c > 0, there exists k independent of u 0 ∈ U 0 such that This follows again by compactness of U 0 in L p , which can be covered by finitely many balls B c/3 (u j ) in L p (T d ). Choosing k large so that we get, for u 0 ∈ B c/3 (u j ), Remark 4.3. One can implement the same scheme to handle more general Hamilton-Jacobi equations of the form where A ∈ C([0, T]; W 2,∞ (T d )) and λI d ≤ A ≤ ΛI d for 0 < λ ≤ Λ. In particular, one has to appropriately adjust the proofs of the integral and Hölder estimates, following [20], and use the linear maximal regularity results in [54].

Applications to Mean Field Games
For a given couple u T , m 0 ∈ C 3 (T d ), consider the MFG system

The monotone (or defocusing) case
Proof of Theorem 1.4. We argue that under the restrictions on r it is possible to prove a priori bounds on second order derivatives of solutions to (MFG) (and beyond, assuming additional regularity of the data). These are typically enough to prove existence theorems. One may indeed set up a fixed-point method, or a regularization procedure, which consists in replacing g(m) by g(m χ ε ) χ ε (where χ ε is a sequence of standard symmetric mollifiers). The existence of a solution (m ε , u ε ) is then standard (see e.g. [31]).
Since bounds on (m ε , u ε ) do not depend on ε > 0, it is therefore possible to pass to the limit and obtain a solution to (MFG).
The key a priori bound is stated in the next Lemma 5.1. Once bounds for u in W 2,1 q (Q T ), q > (d +2)/γ , are established, one can indeed improve the estimates via a rather standard bootstrap procedure involving parabolic regularity for linear equations. Indeed, D p H(x, Du) turns out to be bounded in L p for some p > d + 2, which is the usual Aronson-Serrin condition yielding space-time Hölder continuity of m on the whole cylinder (see e.g. [36, Theorem III.10.1]). We can then use Theorem 1.1 to conclude that u in W 2,1 q (Q T ) for any q > d + 2. This immediately implies by embeddings of W 2,1 q (Q T ) that u is bounded in C 1+δ, 1+δ 2 (Q T ) for any δ ∈ (0, 1). Then, we can regard the Hamilton-Jacobi equation as a heat equation with a space-time Hölder continuous source, and by [36, Section IV.5.1] conclude that u ε is bounded in C 2+δ , 1+δ 2 (Q T ) independently of ε. One then goes back to the Fokker-Planck equation to deduce that m also enjoys C 2+δ , 1+δ 2 (Q T ) bounds.
Below we state and prove the crucial a priori estimate on solutions to (MFG).
Step 2. Second order estimates. These are obtained by testing the Hamilton-Jacobi equation with ∆m and the Fokker-Planck equation with ∆u: On one hand, (H2), while on the other hand by Young and Cauchy-Schwarz inequality therefore, back to (39), integrating by parts we obtain Plugging in (38), using lower bounds on g and the fact that u ≥ min u T + min H(·, 0) by the comparison principle, we finally get Step 3. Setting b(x, t) = −D p H(x, Du(x, t)), by the assumptions on H, (38) and (40)  if γ > 2.
Under our assumptions on r, the exponent q can be chosen large enough to apply Theorem 1.1, that yields the assertion.
The following lemma is needed to extract regularity information on m. Then, there exists a constant C depending on K,µ and d, such that Proof. The estimate can be obtained by testing the equation against m p−1 and using parabolic interpolation (see [30,Theorem 4.1] for further details). We briefly sketch it here for completeness. Testing the equation and integrating by parts yields Then, using also Young's inequality, (41) is less than or equal to . If d > 2, then 2 d(µ−1) < 1 and we are done. If d ≤ 2, the proof is somehow simpler, by different nature of Sobolev embeddings of W 1,2 (T d ); we refer to [32].

The non-monotone (or even focusing) case
Proof of Theorem 1.5. We detail only the derivation of a priori estimates for smooth solutions to (MFG), i.e. the existence of a constant C > 0 (depending only on the data) such that Since we fall in the maximal regularity regime for the Hamilton-Jacobi equation, the existence of a solution to the MFG system can then be derived as in Theorem 1.4. We start from first order estimates as in Lemma 5.1, that now yield, using the assumptions on g, To bound the two terms in the left hand side separately, a further step is needed. Recall a crucial Gagliardo-Nirenberg type inequality proven in [22,Proposition 2.5], that reads |D p H(x, Du(x, t))| γ dxdt + 1 .
for some C > 0 and δ < 1, provided that r < γ /d. Then, using this inequality, the assumptions on L and the fact that Hence, back to (43), we get m L r +1 (Q T ) ≤ C 2 .
Therefore, g(m) L r +1 r (Q T ) ≤ C 2 . Note that by the assumptions on r, q = r+1 r is large enough to apply Theorem 1.1, and (42) follows.
Proof. We just need to localize analogous results on R d , which go back to [61,Lemma 9 p. 441], see also [59, Proposition 10-(b) Section 5.2] and [40,Theorem 17.38] for a proof via real interpolation methods. Let χ be a compactly supported cut-off function such that χ ≡ 1 on the unit cube [−2, 2] d . It is immediate to see that the extension operator is bounded for all nonnegative integers k ≥ 0 and p ≥ 1, so it is bounded from W α, p (T d ) to W α, p (R d ) by interpolation. Then, where the second inequality relies on the aforementioned embedding W α, p (R d ) → N α, p (R d ) B α p∞ (R d ).