Niclas Str6mberg

Computer rmoelastic wear problems

The present paper concerns the numerical treatment of thermoelastic wear problems. The governing equations of thermoelasticity coupled to Signorini contact, Coulomb's friction and Archard's wear are formulated as a system of discrete equations. This equation system is solved, using a Bouligand differentiable Newton method, for five problems of didactic nature.

Introduction

In the present paper, a Newton method for thermoelastic wear problems is suggested. The continuum thermodynamical wear model presented in [START_REF] Strrmberg | Derivation and analysis of a generalized standard model for contact, friction and wear[END_REF] is considered for a two-dimensional thermoelastic body which is unilaterally constrained by a rigid flat tangential moving support. In addition, it is assumed that no heat can flow through the support. For this system, starting from the governing equations, a system of finite element equations, which can be associated with an augmented Lagrangian formulation (see e.g. [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF][START_REF] Klarbring | Mathematical programming and augmented Lagrangian methods for frictional contact problems[END_REF]), is derived. This system of equations is numerically solved for five problems of didactic nature using a Bouligand differentiable Newton method given by Pang [START_REF] Pang | Newton' s method for B-differentiable equations[END_REF].

Augmented Lagrangian formulations of two-and three-dimensional elastic wear problems were solved in [START_REF] Strrmberg | An augmented Lagrangian method for fretting problems[END_REF][START_REF] Strrmberg | A Newton method for three-dimensional fretting problems[END_REF] using Pang's Newton method. Furthermore, an extensive comparison of Pang's algorithm with an interior point method and a method commercially available in ABAQUS was performed for two-dimensional friction problems in [START_REF] Christensen | Formulation and comparison of algorithms for frictional contact problems[END_REF]; numerical tests showed that the Newton method is both effective and robust. In the present paper, thermal properties are added to the augmented Lagrangian formulations, discussed in the papers above, in order to achieve a robust and efficient algorithm for thermoelastic wear problems. The method achieved is a direct method without any splitting between mechanical and thermal subproblems (see e.g. [START_REF] Simo | Associative coupled thermoplasticity at finite strains: formulation, numerical analysis and implementation[END_REF]). Numerical tests have proven that the method works viably. However, a crucial detail for making the method effective is to eliminate the temperature by algebraical manipulations, such that the temperature is not explicitly involved in the Newton step, but rather appears as an inner variable in the algorithm. Otherwise, the temperature field will not be solved accurately, which in turn, of course, would affect the quality of the mechanical part of the solution.

During the past years, a certain degree of progress has been made concerning the finite element analysis of thermomechanical contact problems. ~ Zavarise et al. [START_REF] Zavarise | Real contact mechanisms and finite element formulation--a coupled thermomechanical approach[END_REF] formulated a finite element model for thermal contact conduction; frictional heat generation was briefly discussed using a frictional tangential stiffness term and Concerning the analysis of thermomechanical contact problems of simple geometries such as half-planes, see [START_REF] Barber | Thermoelastic contact problems[END_REF]. 0045-7825/99l$see front matter © 1999 Elsevier Science B.V. All rights reserved. PII: S0045-7825 (98)00392-2 1 additional empirical coefficients. In [START_REF] Wriggers | Thermomechanical contact--a rigorous but simple numerical approach[END_REF], frictionless thermal contact conduction problems were considered using a penalty formulation, where the penalty terms depended on micro-phenomenological assumptions. Thermoelastic contact problems including both thermal contact conduction and frictional heat generation were formulated and solved in [START_REF] Johansson | Thermoelastic frictional contact problems: modelling, finite element approximation and numerical realization[END_REF] by splitting the problem in one mechanical part and one thermal part; the mechanical problem was solved using a frictional contact algorithm suggested by Klarbring and Bjrrkman [START_REF] Klarbring | A mathematical programming approach to contact problems with friction and varying contact surfaces[END_REF]. Wriggers and Miehe [START_REF] Wriggers | Contact constraints within coupled tbermomechanical analysis--A finite element model[END_REF] treated frictional contact problems including thermal contact conduction and frictional heat generation using a normal compliance law and a radial return approach for the friction. In [START_REF] Zavarise | On augmented Lagrangian algorithms for thermomechanical contact problems with friction[END_REF], augmented Lagrangian techniques were used for thermomechanical frictional contact problems described with microscopical laws; unfortunately, as it seems (see also [START_REF] Wriggers | Application of augmented Lagrangian techniques for non-linear constitutive laws in contact interfaces[END_REF]), the theoretical contact formulation is for bilateral contact and not for the proper unilateral formulation; in addition, a simplification of Coulomb's law of friction is used and the frictional heat generation is excluded. Recently, finite element formulations of thermomechanical frictional sliding in large deformations were given by Laursen and Oancea [START_REF] Laursen | On the constitutive modeling and finite element computation of rate-dependent frictional sliding in large deformations[END_REF][START_REF] Oancea | A finite element formulation of thermomechanical rate-dependent frictional sliding[END_REF].

In this paper, we focus on the finite element treatment of frictional heat generation. Therefore, thermal contact conduction is not considered; of course, it is straightforward to include heat transfer through the contact interface (see [START_REF] Strrmberg | Derivation and analysis of a generalized standard model for contact, friction and wear[END_REF]). Actually, wear heat generation is also included in the thermomechanical model presented herein. However, for physically realistic wear coefficients, the wear heat generation can be neglected in comparison to the frictional heat generation. Still, the internal state variable measuring the wear gap has a major influence on the frictional heat generation and vice versa. This fact is illustrated by numerical examples. For instance, it is numerically shown that two frictional hot spots appearing in a frictional analysis are removed if wear also is included in the constitutive settings.

The contents of the present study is as follows: in Section 2, the governing equations for a unilaterally constrained thermoelastic body are given. In Section 3, the governing equations are put together to an initial boundary value problem, which in turn is discretized both in space and time in order to reach a system of finite element equations. The system of equations obtained is considered in Section 4, where the adoption of Pang's Newton algorithm also is discussed. In Section 5, five numerical problems are studied. The problems differ only in the constitutive settings. That is, five different combinations of elasticity, thermoelasticity, Signorini contact, Coulomb's friction and Archard's wear are investigated for a particular initial boundary value problem. Finally, in Section 6, concluding remarks are given.

Governing equations

Let us consider a two-dimensional body in a quasi-static small displacement setting, occupying a domain /2 @ ~t2 with a piecewise smooth boundary a~O, see Fig. the heat flux vector by q, the Cauchy stress by 0", the infinitesimal strain tensor by ~ and the location of a point by x. Vectors and tensors are represented by components given in a Cartesian basis (e l, e2) which is defined in Fig. 1. Moreover, the summation convention is applied to repeated indices and Kronecker's delta 6,j is used. The boundary 0/2 is dived in three disjoint parts/7, F, and ~. The body is subjected to prescribed tractions t on /7, and zero displacements on F . It is further assumed that the body is unilaterally constrained on /7,. by a rigid support in the direction n,.. No heat transfer can take place through the support. Furthermore, the support is flat ended with a normal n s = e 2 and subjected to a tangential displacement ~ = sCe~, where s c = ((t) is a function of time t. An important consequence of the small displacement assumption is that n, = -n , . This fact is utilized in the work by setting n, = -e 2.

Equilibrium equations and principles of thermodynamics

Neglecting inertial forces and body forces, the equilibrium equations for the system considered above read ao; v.

= 0

( i = 1 , 2 ) i n / 2 , (1) axj o'ijn j = t i ( i = 1 , 2 ) on/7, (2)
°12 = P t , 0"22 = -P , on ~, (3)
where n i is the i:th component of the outward unit normal on /7, and p =p,e~ -p , e 2 represents the contact traction vector acting on the support. Two work conjugate measures belong to p. These are u, = -u 2 and w, = u~ -s ~. That is, u n is the normal displacement and w, is the relative tangential displacement between the contact surface and the support.

The balance of energy of the body and the contact surface are given by

Oqi pd = o'g/~i/ Oxi in a2, (4)
/~ = p ~/ + p , ~b , -qzon ~, [START_REF] Ghasemi | Surface temperatures and fretting corrosion of steel under conditions of fretting contact[END_REF] where p is the mass density of the bulk material, e is the specific energy of the body, E is the density of the internal energy of the contact surface and a superimposed dot represents a time derivative. Defining the absolute temperature T, the specific entropy of the body s and the entropy density of the contact surface S, the second principle of thermodynamics is given by the following Clausius-Duhem inequalities:

p ~ <~ 0-ij~i; -psi" -q~ aT T Ox i in g2, (6)
~t <~p, ti, + p,ffe -SJ" on ~,

where gt = e -sT and q* = E -ST are the Helmholtz free energies of the body and the contact interface, respectively.

Constitutive assumptions--body

Established constitutive assumptions of isotropic thermoelasticity are utilized for the body. 2 That is, defining Lame's elasticity coefficients .~ and /2, the specific heat capacity c, the thermal dilatation coefficient a and a reference temperature T n, the following free energy is considered for the body:

1 - 2 c = ~p (a(g,) + 2/2~je~i -(6.~ + 4/2)cr(T -To)g/) -~-~o(T -To) 2, (8)
from which the following state laws are derived per definition:

z The material given in this section is classical and can be found in textbooks (see e.g. [START_REF] Lemaitre | Mechanics of Solid Materials[END_REF][START_REF] Maugin | The Thermomechanics of Plasticity and Fracture[END_REF]).

In agreement with this latter inequality, representing the second law of thermodynamics, Fourier's law of heat diffusion is adopted. That is,

OT qi=-k,~x i (i=1,2), ( 13
)
where k, is the thermal conductivity. Inserting (2.13) in the balance of energy represented by [START_REF] Klarbring | Mathematical programming and augmented Lagrangian methods for frictional contact problems[END_REF] yields

k, ~x~ = (3~ + 2~z )aT~kk +-~o TT . ( 14
)
This equation is rewritten, by neglecting ~kk and putting T = T o, to read k, Tx, =ocT, [START_REF] Oancea | A finite element formulation of thermomechanical rate-dependent frictional sliding[END_REF] which is the standard heat conduction equation of uncoupled thermoelasticity.

Constitutive assumptions--contact surface

The constitutive assumptions of the contact surface are given by Signorini's contact conditions, Coulomb's law of friction and Archard's law of wear. These laws were derived from a generalized standard material, by defining a certain internal state variable o~ measuring the wear gap, in [START_REF] Strrmberg | Derivation and analysis of a generalized standard model for contact, friction and wear[END_REF]. Utilizing this internal state variable approach, the free energy is assumed to be given by qt=Ic(u,,,o9 ), C={(u,,og):u,-oo-g<~O}, [START_REF] Pang | Newton' s method for B-differentiable equations[END_REF] where I c is the indicator function of the closed convex set C (i.e.

I c = 0 if (u,,w) E C, 1 c = +oc if (u,,, o~) ~ C)
and g is the initial gap between the body and the support. In addition, the following state laws are defined and derived:

=fNc(u., w) if (u.,~o) EC (p,,, -7~) ~ O!g(u,, O3) V otherwise, ( 17
) 0!/' -s = -b--~-= o, ( 18
)
where °/4/" is the force associated to o~, called the wear driving force. O~ stands for the subdifferential of ~ and

Nc(u ,, w) = {(p,, -7/t/') : 0 >~ p,,(u~ -u,,) -7/U(w' -w) V (u~, w') E C}

defines the normal cone of C. The contact law given in [START_REF] Petiot | An analysis of fretting-fatigue failure combined with numerical calculations to predict crack nucleation[END_REF] represents an extension of Signorini's contact conditions that includes wear.

To evaluate the implication of the choice of free energy and state laws given in ( 16)-( 18) in conjunction with the balance of energy in [START_REF] Ghasemi | Surface temperatures and fretting corrosion of steel under conditions of fretting contact[END_REF] and the Clausius-Duhem inequality in [START_REF] Hughes | The Finite Element Method, Linear Static and Dynamic Finite Element Analysis[END_REF], we consider evolutions such that the left-hand time derivative and the right-hand time derivative of the state variables and the state functions exist with possibly different values, meaning that the time derivative of these functions and variables may not necessarily exist. That is, the state functions and the state variables are continuous in time, but not necessarily smooth. Thus, a superimposed dot stands for either the left-hand time derivative or the right-hand time derivative.

Then, the above assumption on the free energy implies that ~= 0 for (u., w)~ C; otherwise ~ is not defined. Besides, (u., w) ~ C implies that ~ = +~ which indicates that such states are impossible. Moreover, [START_REF] Simo | Associative coupled thermoplasticity at finite strains: formulation, numerical analysis and implementation[END_REF] implies that S = 0. Consequently,/~ = ~ + $7 ~ + ST = 0 for thermodynamical admissible processes, i.e. for (u., w)E C. Furthermore, the state law given in ( 17) is equivalent to 7g'=p,>~O, u~-w-g<-O, p,(u,-w-g)=O.

(

) 19
These conditions imply also that p~i =p,&. Putting ~=/~= S = 0 and p,d, =p,& in ( 5) and ( 7) yield

q2 = p d~ + p,vb,, (20)
p, o3 + p,#,/> 0.

These two relationships represent the balance of energy and the second principle of thermodynamics for the contact interface. Furthermore, for later use, we remark that the contact conditions in [START_REF] Sproles | Interface temperature measurements in the fretting of a medium carbon steel[END_REF] are optimality conditions to the following variational inequality:

p.>~O: (u~-w-g)(p~-p,,)<-O Vp'~>0. ( 22
)
Next, Coulomb's law of friction and Archard's law of wear are introduced. Coulomb's law of friction is given by the following principle of maximal dissipation:

IP, I~</~(P,,)+ : G(P~ -p,)<-O V[pil~</z(p,,)+, (23)
where/z is a constant friction coefficient, I1 stands for the absolute value and (p.)+ = (p. + Ip.[)/2. Assuming that k is a constant wear coefficient, Archard's law of wear is represented by

;, = k(p.)+ Iw,I, (24)
One might remark that the use of (p.)+ in ( 23) and [START_REF] Waterhouse | Fretting Corrosion[END_REF], instead of just p., implies that these two representations of Coulomb's and Archard's laws are defined also for p. < 0 and not only for p,,/> 0.

Let us evaluate the consequences of Coulomb's law and Archard's law with respect to the principles of thermodynamics given in [START_REF] Strrmberg | Derivation and analysis of a generalized standard model for contact, friction and wear[END_REF] and [START_REF] Strrmberg | An augmented Lagrangian method for fretting problems[END_REF]. Taking p~ = 0 in [START_REF] Vingsbo | On fretting maps[END_REF] implies that ptvb, >/0. Using this and [START_REF] Waterhouse | Fretting Corrosion[END_REF], it is shown that the second law of thermodynamics represented by ( 21) is satisfied. Furthermore, [START_REF] Vingsbo | On fretting maps[END_REF] implies that p, =/z(p,)+sgn(vb,) when w, ¢ 0. This together with [START_REF] Waterhouse | Fretting Corrosion[END_REF] in [START_REF] Strrmberg | Derivation and analysis of a generalized standard model for contact, friction and wear[END_REF] result in the following energy balance: q2 :-(k(P,,)2+ +/x(P,,)+)lwtl, [START_REF] Wriggers | Contact constraints within coupled tbermomechanical analysis--A finite element model[END_REF] which represents the friction and wear heat generation at the contact interface. Thus, all heat generated by friction and wear flows into the body.

The system of discrete equations

Given proper initial conditions, a load history t(t) and a displacement history ~(t) on a time interval [0, r], the following initial boundary value problem is stated: Find u :[0, z] ~ 7/', p. :[0, ~']----) K,,, p, :[0, z]---)F(p.) and T : [0, 7] --4 3-such that for each t E [0, ~-]: and 3-is a function space of temperatures. In addition, weak forms of the tribological laws given in ( 22)-( 24) are assumed to be valid on F c (see [START_REF] Strrmberg | An augmented Lagrangian method for fretting problems[END_REF]). Eq. ( 26) is the principle of virtual work, equivalent to the equilibrium equations given in ( 1)-( 3), with the thermoelastic properties in (9) inserted. Finally, Eq. ( 27) is the balances of energy for the body and the contact surface given in ( 15) and ( 25) put together using [START_REF] Lemaitre | Mechanics of Solid Materials[END_REF] and q.n = -q2.

f. ov, oo; fr Ely ax~-~ O~ dV-G(T -To) ~x~ dV=

tiv; dA -, (p,v, -p.v2) dA V v E 7#, ( 26
)
pcTOdV+ k,-~x~xdV= ( (p,,)+ +~(p,,)+)lw,lOdA roe3-, (27)

Space discretization

The initial boundary value problem stated above is treated by introducing finite element approximations of 7/" and if (see e.g. [START_REF] Hughes | The Finite Element Method, Linear Static and Dynamic Finite Element Analysis[END_REF]), and evaluating the integrals over I~ by an appropriate quadrature rule depending on the choice of finite elements. Let f be a function over F,, then the integrals are approximated by

fry da ~-- IMf(xM), ( 28
) Z METI¢.
where I M are the weighting factors and r L is the set of integration points x M on 4-When applying this formula, the contact boundary is divided into contact elements. If these elements coincide with the edges of the displacement elements, then it is possible to choose an integration rule such that the integration points will coincide with the nodal displacement points of the contact surface. For two-dimensional problems, examples of such rules are the trapezoidal quadrature rule for linear elements and Simpson's rule for quadratic elements. A study of different integration rules can be found in [START_REF] Johansson | A study on the discretization of a mixed formulation of the frictionless contact problem[END_REF].

After performing the space discretization outlined above, one obtains the following discrete equilibrium equation from [START_REF] Wriggers | Thermomechanical contact--a rigorous but simple numerical approach[END_REF]:

/i'd -/~(T -T o) = F T -C,,P. -CV, P, ( 29
)
where d = {d A} represents the nodal displacements, T = {T A} is the nodal temperatures and T o = {To A} is the nodal reference temperatures at the nodal points x A. Furthermore, letting N A= NA(x) represents the shape functions,

-BA -ea fo oNA ON----~e g =[Kik ]' Kik =_..Eij~l Ox I Oxj dV, 1~ ^ BA ^ eA fa ONe = [K i ], K i GN a ~ dV, F = {F~} FB = Jr, tiNB dA ' i MB MB C, = [C,i ], C,, i = --Ne(xM)62i , MB Me B M C, = [C,i ], C,i = N (x )~1i , it',, {P,7} pM M M = , ,=Ip,(x ), p, . M = , =I p,(x ),
The discretization of the balance of energy in [START_REF] Wriggers | Application of augmented Lagrangian techniques for non-linear constitutive laws in contact interfaces[END_REF] can be expressed by

MT + OT = L(a, P,), (30)
• M k to = ~ (P~)+l vb,~l.

(33)

Here, u~ = --NA(xM)d A, toM = to(X M) and gM = g(xM). In [START_REF] Strrmberg | An augmented Lagrangian method for fretting problems[END_REF], the above relationships were derived from weak forms of ( 22)-( 24) by using the integration rule in [START_REF] Zavarise | On augmented Lagrangian algorithms for thermomechanical contact problems with friction[END_REF]. Finally, we rewrite (31) in the following manner for any r > 0 (see [liD: find P,~ t>0 such that

(p,M _ (py + r(u y _ tom M M' M M' -g )))(P. -P , . ) ~O VP,, 9 0 , ( 34
)
implying that pM is the projection of pM + r(u,M --toM --gM) onto K h,, = {p,M: p,M/> 0}.

Time discretization

The time rates appearing in (30), (32) and (33) are approximated using a backward Euler discretization, i.e.

T(tk + 1 ) --T(tk + l ) -T(tk)

tk+ i --tk ' (35)
(M + (tk+ l -tk)O)T(tk+l) = Q + £(d(tk+l).P (tk+j)). ( 38
)
where Q = MT(tk) and

f ~B M / k M 2 ) } L ( d . P . ) = i & N ( x ) ~-f f ( P )+ +#(PM)+ i ~1 . ( 39
)
By using (36) and the definition of ~M. Coulomb's law in (32) can be rewritten, for any r > 0 . as the M M following variational principle: find ]P, I <~/z(P,. )+ such that -M M '

(pM _ (pM + rw, ))(Pt -pM) >/0, V Ie7'l <~/x(P,M)+, where to,M = toM(tk)"

(41)

The algorithm

The initial boundary value problem defined in the beginning of Section 3 can now be treated by solving a sequence of the discrete equations given in ( 29), (34), (38), ( 40) and (41) using Pang's algorithm presented below. However, numerically tests have proven that the search direction involved in the algorithm can not be solved accurately for this representation of the equations. Comparisons by analytical solutions given in [START_REF] Carslaw | Conduction of Heat in Solids[END_REF] have shown that the Newton step smoothens out the temperature gradients in a wrong way; in addition, the algorithm converges slowly. This depends probably on the great differences in the numerical values of /~, /~ and (M + (tk+ j -t~)O), leading to a system of ill-condition equations. This problem is overcome by putting together [START_REF] Zavarise | Real contact mechanisms and finite element formulation--a coupled thermomechanical approach[END_REF] e,,TP,,e,, (42) where /~ =/~(M + (tk+ j -tk)O) 1. In such manner, the temperature is eliminated from the Newton step and appears instead as an inner variable in the algorithm.

and (38) to form

Thus, letting y = (d, P,,, P,), the following system of equations are solved for each time increment:

- _~_ T -t-T )

Kd-lf(Q+L(d,P,,)) c,,e, C,P,+KTo-F H(y) = -P,, + H,,(d, P,,) -P, + ll,(d, P,,, e,)

= 0, ( Fh.p Mis the projection of P~ + r~M onto ( , ). The equation system in ( 43) is an extension of the quasi-augmented Lagrangian formulations considered previously in [START_REF] Strrmberg | An augmented Lagrangian method for fretting problems[END_REF][START_REF] Strrmberg | A Newton method for three-dimensional fretting problems[END_REF]. These formulations are ~imply obtained by setting /(" = 0 in the above formulation.

Pang's Newton algorithm

In [START_REF] Strrmberg | An augmented Lagrangian method for fretting problems[END_REF], it was shown that, for/~ = 0, the equation system in (43) is Bouligand differentiable, i.e. Lipschitz continuous and directionally differentiable. Obviously, this holds also for the case with /~ # 0 as/~(d, P,,) is Lipschitz continuous and directionally differentiable. For solving systems of Bouligand differentiable equations, Pang [START_REF] Pang | Newton' s method for B-differentiable equations[END_REF] suggested the following globally convergent algorithm (concerning convergence result of the algorithm, see also [START_REF] Christensen | Formulation and comparison of algorithms for frictional contact problems[END_REF]):

ALGORITHM. Let fl, y and • be given scalars with fl E (0, 1), y E (0, 1/2) and e small. Repeat the following steps for each time increment k + l" (0) Let yO be given from the previous time step k and set l = 0.

(1) Find a search direction z = (za, z,,,zt) such that

H(f) + H'(f; z) = 0, ( 46
) t I
where H (y ;z) is the directional derivative. (2) Let a ~= tim,, where m~ is the first non-negative integer m for which the following decrease criterion holds:

g(f + fl~Z) <~ (1 -2yflm)g(f),

T g(y) = ~ H (y)H(y).

(3) Sety I+l =f + dz.

(4) If g(f+l) ~< e, then terminate with y/+ 1 as an approximate zero of H(y). Otherwise, replace l by 1 + 1 and return to step 1. After convergence in each time step is achieved, the variables T, w~ and w M are updated.

Approximation of directional derivatives

The time-consuming part of the algorithm is to find the search direction in (46). Generally, this system of equations is non-linear due to the non-linearity of H'(y; z) in the second argument, which in turn depends on the non-differentiability of H(y). Of course, H'(y;z)=VH(y)z at states where H(y) is differentiable. The nonlinearity of (46) might limit the effectiveness of the algorithm. This drawback is overcome by introducing an appropriate approximation of H'(y;z) such that (46) results in a system of linear equations (see [START_REF] Christensen | Formulation and comparison of algorithms for frictional contact problems[END_REF][START_REF] Strrmberg | An augmented Lagrangian method for fretting problems[END_REF][START_REF] Strrmberg | A Newton method for three-dimensional fretting problems[END_REF]).

The non-differentiability of H(y) concerns the functions /2(y)=£(d,P,,), /~,,(y)=//~,,(d,P,,) and //~,(y)= /~,(d, Pn, P,); otherwise, H'(y; z)= VH(y)z. In order to obtain a linear equation system in (46), the following approximations of/2'(y; z), lI~(y; z) and II[ (y; z) are introduced: :P,,[w,I<~O,P,,+r u,,,,)+[~,[ :&lw, l>o,& +r u,,-,Oo -~e,,l~,l-g ~ >0},

f x-, 8 M I/{2kP,M, ) {kp,M, )~M-M M\a £'(Y;Z)=~ 2., N (x)l|---~-+/z [w,

Numerical examples

Five different problems are considered for a structure subjected to boundary conditions and a displacement history as given in Fig. 2. The problems differ only in the constitutive settings, i.e. the same boundary conditions and displacement history (except for the number of cycles) are used for all five problems. The dimensions of the structure are 3.0 X 1.0 cm 2. Furthermore, the structure is fixed at the top and unilaterally constrained at the bottom by a rigid support using a negative initial gap of g = -0.1 txm. In addition, the support is subjected to a cyclic tangential displacement history with an amplitude of 1.0 mm according to Fig. 2. In the context of wear, this type of boundary condition can be called a fretting condition. Although ~: is small, this condition will result in global sliding fretting. It is a major principle difference between the contact stresses obtained for global sliding fretting compared to stick-slip fretting. This is discussed in more detail further on. The thermoelastic properties of the structure are approximated by using 42 X 28 bilinear quadrilateral elements. [START_REF] Carslaw | Conduction of Heat in Solids[END_REF] The potential contact surface at the bottom of the structure is approximated by 42 contact elements using the trapezoidal rule in such manner that the integration points coincide with the nodal displacement points of the contact surface. Concerning the time discretization, each cycle is treated using four time steps, i.e. tk+ j = t k + 0.01 s.

For the system given above, the following problems, defined by five different constitutive settings, are considered: the first problem is pure elastic, meaning that no thermal properties are used, with /.t = k = 0, the second problem is pure elastic with /~ = 0.3 and k = 0, the third problem is pure elastic with /x = 0.3 and k= 1.0x 10 -tj Pa -~ (this is a typical wear coefficient for adhesive wear in steel assemblies), the fourth problem is thermoelastic with # = 0.3 and k = 0, and, finally, the fifth problem is thermoelastic with /z = 0.3 and k = 1.0 X 10-~ Pa -j.

The contact pressures of problem one (/1, = k = 0, no heat) and problem two (/x = 0.3, k = 0, no heat) are compared in Fig. 3. The only difference between these two problems is the value of the coefficient of friction. Obviously, p,, is constant in time for/z = 0, but, for # = 0.3, p,, depends on the sliding direction of the support. In Fig. 3, p,, is plotted at states when t = n X 0.04 s, where n is a positive integer. The solution will be the same at times when ~ > 0 and global sliding is developed. On the other hand, the reflected image of p,, will appear when there is global sliding with ~ < 0. Furthermore, the dissipation is non-zero for/z = 0.3 in contrast to the case for /z = 0 when it is zero. Actually, the frictional dissipation is constant for each cycle. Thus, if one is interested in studying the dissipation for a large number of cycles such as the dissipation in frictional dampers, then the constitutive assumptions of problem two might not be sufficient due to the fact that real engineering materials eventually will be broken during a repeated frictional process and consequently the dissipation will eventually be equal to zero. This fact is included in the constitutive settings of problem three and five. For those problems the frictional dissipation turns toward zero after many cycles depending on the surface break down due to wear.

In Fig. 4, the contact pressure and the wear gap are plotted for problem three (/z = 0.3, k = 1.0 X 10-j~ Pa -~, no heat) after 30 cycles (t = 1.2 s). At this instant, the contact pressure is approximately zero since the initial gap is almost worn away, i.e. w ~ -g = 0.1 ;xm over the contact surface. However, during the first cycle of this fretting process, the wear gap is approximately zero and the contact pressure is similar to the pressure obtained in problem two, see Fig. 3. Thus, although a physically realistic value on the wear coefficient is used, a large reduction of the contact pressure has been obtained within a few number of cycles. This depends not only on the global sliding conditions but also on the fact that the wear model used in this work does not include the transportation of wear debris in the contact interface but rather annihilates the wear debris. Of course, the rate of reduction in contact pressure can be decreased by diminishing the wear coefficient. Perhaps a better approach to obtain such a reduction would be to extend the formulation used in this work by a third body of wear debris, following the ideas of Godet [START_REF] Godet | The third body approach: a mechanical view of wear[END_REF]. Another approach in that direction is accomplished in problem five, where the rate of reduction in contact pressure is decreased by adding thermal properties to the constitutive settings of this problem.

One might notice that the result obtained above would be completely different if the fretting conditions would be chosen such that no global sliding would be developed. If the amplitude on ~: would be taken sufficiently small such that no slip would be developed over the contact surface (global stick), then, of course, no wear would be initiated. However, in the intermediate with amplitudes on ( implying stick-slip, the contact pressure would turn toward zero in the slipping regions and the pressure would increase in the sticking regions (see [START_REF] Strrmberg | An augmented Lagrangian method for fretting problems[END_REF][START_REF] Strrmberg | A Newton method for three-dimensional fretting problems[END_REF]). In addition, large values and large gradients in the contact stresses would be developed in the regions between stick and slip. Concerning the field of fretting fatigue, such stress distributions might initiate cracks near the contact surface leading to fatigue failures of the structure. Similar observations have been found experimentally by using so-called fretting maps in e.g. [START_REF] Petiot | An analysis of fretting-fatigue failure combined with numerical calculations to predict crack nucleation[END_REF][START_REF] Vingsbo | On fretting maps[END_REF].

x 10 -5 11=0.3, k=10 -ll, no heat

x 10 -7 la--0. 11

In Fig. 5, two frictional hot spots and two peaks in the contact pressure are plotted after 30 cycles for problem four (/z = 0.3, k = 0, heat). After the first cycle the contact pressure is similar in shape to the pressure of problem two, see Fig. 3, but due to the thermal expansion induced by frictional heat generation the contact pressure changes rapidly from this shape during the evolution of cycles. The shape obtained in the temperature field and the contact pressure are very sensitive to the geometry and the boundary conditions. For instance, if the dimensions of the structure are changed to 1.0 × 1.0 cm 2, then only one frictional hot spot will appear. Furthermore, if the velocity of the support is decreased sufficiently, then no hot spots at all will appear. An obvious drawback of this constitutive model, concerning this type of evolution problems, is that no damages are induced due to the high temperatures and high contact stresses, neither in the bulk material nor in the contact surface. In real engineering materials such high temperatures and high stresses would result in damages of the material in some way, which in turn would reduce the amount of frictional heat generation and consequently result in lower temperatures and pressures than predicted here. This fact is included in the constitutive settings of problem five by adding Archard's law of wear to the constitutive laws used here. In that problem, the surface damages in form of wear imply that no frictional hot spots are developed.

In Fig. 6, the temperature field, the contact pressure and the wear gap are plotted after 210 cycles (t = 8.4 s) for problem five (/x = 0.3, k = 1.0 × 10 ~ Pa -~, heat). At this instant, the solution has almost reached a steady state. The temperature field is approximately constant at 0.0987 K. Notice, in comparison to problem four, that no frictional hot spots are developed. Instead the increase in temperature is astonishingly low. This fact depends on the wear process, which decreases the frictional and wear heat generation. However, following the discussion on problem three, this reduction is probably too large due to the annihilation of wear debris. If the wear debris would be included in the model as a third body, then the frictional and wear heat generation would be higher, but not as high as in problem four. Although the temperature is very low, it has a major effect on the wear rate. After seven times more cycles compared to problem three, the contact pressure is still much higher here. This depends on the thermal expansion. Furthermore, the thermal expansion also increases the amount of wear; the wear gap is no longer constant and the maximum in wear gap obtained is more than 20% larger in comparison to problem three. The temperature rise in fretting is a subject of considerable interest. Some researchers claim that the presence of et-Fe20 3 in the wear debris must be taken as evidence that the interface temperature exceeds 300-600 K in fretting (see e.g. [START_REF] Waterhouse | Fretting Corrosion[END_REF]). Obviously, such high temperatures are not predicted in problem five. Nevertheless, by neglecting the removal of material, as in problem four, even higher temperatures are predicted, see Fig. 5. However, measurements by other researchers indicate that the average temperature in contacting asperities subjected to fretting conditions is much lower; magnitudes of 20-30 K have been measured (see e.g. [START_REF] Ghasemi | Surface temperatures and fretting corrosion of steel under conditions of fretting contact[END_REF][START_REF] Sproles | Interface temperature measurements in the fretting of a medium carbon steel[END_REF]). Such temperature increases are more comparable to the increases obtained by the wear model used in the present paper.

For fretting problems with stick-slip conditions, the frictional heat generation has a minor influence on the contact pressure and the wear gap using this wear model. Still, other initial and boundary conditions on the temperature and the heat flux vector can have a major effect on the contact pressure and the wear gap for stick-slip fretting problems. For instance, one might think of problems where external thermal boundary conditions induce stick-slip fretting. This type of problem can of course be treated with the algorithm suggested in this work.

Finally, some execution statistics for the five problems considered are provided. The number of Newton iterations per time increment (New./inc.) and the number of line-searches, i.e. Armijo iterations, per Newton iterations (Arm./New.) are given in Table 1. In order to get a fair comparison, the statistics are determined for the ten first cycles of each problem. Apparently, problem four is the most time-consuming to solve, and problem one and two are the fastest ones. The reason for problem two to have the same statistics as problem one is that there is only global sliding and no stick-slip developed in the solution of this problem as well as in the solutions of the other problems. If there would be stick-slip in the solutions, then the statistics would be higher for problem 2-5, at least the number of New./inc. would be higher.

Concluding remarks

In this work thermoelastic wear problems have been solved by a finite element method. The method is a direct method, meaning that no iterative strategy between a mechanical and a thermal part of the problem is used, but instead the finite element equations are solved using a Bouligand differentiable Newton method. The method obtained is effective and robust.

Numerically, it has been found that the internal state variable measuring the wear gap has a great influence on the temperature field generated by frictional heat. For a particular frictional heat generation problem, it is numerically shown that when this internal state variable is included in the analysis, then the increase in temperature is very low. On the contrary, if wear is excluded, then two frictional hot spots are developed for the same problem. However, although the temperature increase is very low for the wear problem, the frictional heat generation has a major effect on the wear rate.

The wear model utilized in this work does not include the transportation of wear debris, but instead the internal state variable measuring the wear gap rather annihilates the wear debris. From the numerical result discussed above, it is clear that the transportation of wear debris will affect the frictional heat generation which in turn will affect the wear rate. Therefore, for future plans, it would be of interest to extend this wear model by a third body of wear debris and to study its influence on the wear rate.

1 .Fig. 1 .

 11 Fig. 1. The two-dimensional body considered.

  -~%. = ,~Ekk~,~ + 2/iEsj --(3~ + 2/i)oz(T --To)6~j (i, j = 1, By using d = 0 + s] ? + s'T, and (8)-(10) in (2.4) and (2.6), we obtain Oqi pc 0 = ~ + (3.~ + 2~)aT~kk + -~o TT",

c where 5

 5 E~j,, = ASi/Sk, + 2flzS, k~,, G = (3,~ + 2/2)ce, ~= {v : v(x) = O, x ~ F} , g,, ={p,, : p,,(x) ~> 0, x ~}, F(p,) = {p, : Ip,(x)l <~ #(p,(x)) +, x E 4}

  , P , , ) : { L " } , L " = E N (x )~.~(py)2+ +#( ,,)+ [w,i ME~ e and, in addition, w M = NA(xM)d A --~.Furthermore, the discretization of the tribological laws can be obtained by assuming that (22)-(24) hold in every integration point M E r/,, i.e. I~<~(P.)+ : w, (P, -P ~) < ~o v[PM'I~</x(P,,)+,

  W, (tk+ I) --W, (tk). and inserting (35) and (36) in (30) yield

  Mwhich shows that pM can be expressed as the projection of pM + rv~M onto F (P,,) = {pM: [P,~l ~</x(P, )+} (see[ll]).Finally, using (36) and (37), Archard's law in (33) is approximated by 7 N. Strgmberg / Comput. Methods Appl, Mech. Engrg. 177 (1999) 441-455 Mt k M to (k+,)--~ too M +Tff(P,, (tk+l))+[w~l,

/

  ~,,(d,P~) = pM+r U, -w,, --~(P,,)+lw, [_gM n + (44) is the projection of P~ + r(u, M, -w M inserted and _ gM) onto Kh, with the approximation of Archard's wear law given in (41) , I~/z(P )+'~ II,(d, P,,, P,) = [#(p,,M)+ sgn(P,M+ rw,-M) otherwise' J (45)

  ,)÷]w~I--gM ~<0}, rI,.2 = {M : P, M, r u,, -to,

  , >O,IP, +r~'~]~</-tP,,}, M : P,, >0, ]P7 + r*Yl >/xe,~} 9 and N. Str6mberg / Comput. Methods Appl. Mech. Engrg. 177 (1999) 441-45.5 Zd,, = C,,Zd : {Zd,,} , Zd, = C,Zd : {Zd~} , Z,, : {Z,~} and Z, : {Z •

Fig. 2 .

 2 Fig.2. The structure with the boundary conditions and the displacement history for the five problems considered.

Fig. 3 .

 3 Fig. 3. Comparison of the contact pressures for problem one (tz =k=0, no heat) and problem two (/z =0.3, k=O, no heat).

Fig. 4 .

 4 Fig. 4. The contact pressure and the wear gap for problem three (/x=0.3, k = 1.0× 10 " Pa ', no heat) plotted after 30 cycles.

Fig. 5 .

 5 Fig. 5. Two frictional hot spots and two peaks in the contact pressure are developed after 30 cycles for problem four (tz =0.3, k=0, heat).

Fig. 6 .

 6 Fig.6. The temperature field, the contact pressure and the wear gap plotted after 210 cycles for problem five (tz =0.3, k = 1.0× 10 ~ Pa -~, heat). Notice that the scale on the temperature is given in milli-Kelvin.

13

 13

  ........... ) ............. ) ............ ......................... i. ............

						ILt=0.3, k=10 -ll, heat
	98.77, 98.765~ ' .....'"	..?"	•	• ... ..... ~ ) i .. ..... i	0.3 ~ 0.225 o.~ 0.15
	98.76,					0.075
	~ 98.755, [-, 98.75. 98.745,					0 1.22	0 x 10 -7	0.01 g--0.3, k=lO -ll, heat 0.02 x I [ml	0.03
	98.74, • ......	"" ................. i".. i	1.195
	....-'"		"'..	.....'"	""..... !	'~ 1.17
		•	'"'"::~"..	0.01
	0.02 ~/,I"'0.005
	x l[m] 0.03 0		x 2[m]

i.t__0.3, k=lO-ll, heat .... "'!.. .

Table 1

 1 Some execution statistics for the five problems considered

	Problem-constitutive setting	New./inc.	Arm./New.
	(1) /x =k=0, no heat	1.0	1.0
	(2) /x = 0.3, k = 0, no heat	1.0	1.0
	(3) /x=0.3, k=l.0Xl0 ~Pa -I,noheat	2.25	5.16
	(4) /~ = 0.3, k = 0, heat	2.7	3.96
	(5) #=0.3, k= 1.05< 10-~Pa -t, heat	1.725	3.07

Acknowledgments

This work was supported by the Swedish Research Council for Engineering Sciences (TFR). I wish to thank Anders Klarbring and Lars Johansson for fruitful discussions.


P Jameson

Alp Ár

R M Ész Áros

SOBOLEV REGULARITY FOR FIRST ORDER MEAN FIELD GAMES

Keywords: mean field games, Hamilton-Jacobi equations, Sobolev regularity of the solutions 2010 AMS Subject Classification: 49K20, 35Q91, 49N60, 49N15, 49N70. 1

come

Introduction

The theory of Mean Field Games (briefly MFG in the sequel) saw the light more than a decade ago, thanks to the works of J.-M. Lasry and P.-L. Lions on the one hand (see [START_REF] Lasry | Jeux à champ moyen I. Le cas stationnaire[END_REF][START_REF] Lasry | Jeux à champ moyen II. Horizon fini et contrôle optimal[END_REF][START_REF] Lasry | Mean field games[END_REF]) and M. Huang, R. Malhamé and P. Caines, on the other hand (see [START_REF] Huang | Large population stochastic dynamic games: closed-loop McKean-Vlasov systems and the Nash certainty equivalence principle[END_REF]). Their main motivation was to study limits of Nash equilibria of (stochastic or deterministic) differential games when the number of players tends to infinity. Since then, it became a very lively and active branch of the theory nonlinear partial differential equations.

In its most simple form an MFG can be described in an informal way as follows. As data, one considers the given quantities: T > 0 (the time horizon), H : T d × R d → R (the Hamiltonian), φ T : T d → R (the final cost of the agents), f : T d × P(T d ) → R (the running cost of the agents) and m 0 ∈ P(T d ) (the initial distribution of the agents), where T d := R d /Z d denotes the d-dimensional flat torus and P(T d ) stands for the set of nonnegative Borel probability measures on T d . A typical agent predicts the evolution of the agents' density m : [0, T ] × T d → [0, +∞) and to find her/his optimal strategy, solves the control problem (1.1) inf α T t L(γ(s), α(s)) + f (γ(s), m(s, γ(s))) ds + φ T (γ(T )) =: φ(t, x)

subject to γ ′ (s) = α(s), s ∈ (t, T ] γ(t) = x ∈ T d .
Here the Lagrangian L : T d × R d → R is the Legendre-Fenchel transform of H w.r.t. the second variable. The optimal strategy is given in feedback form, hence for the agent it is optimal to play -D ξ H(γ(s), ∇φ(s, γ(s))). Having this velocity field as drift, the evolution of the agents' density is given by the solution of a continuity equation. If the prediction coincides with this evolution, the game has a Nash equilibrium. Thus, searching for Nash equilibria for MFG is equivalent to solving the following system of PDEs:

(1.2)

     -∂ t φ + H(x, ∇φ) = f (x, m), in (0, T ) × T d , ∂ t m -∇ • (D ξ H(x, ∇φ)m) = 0, in (0, T ) × T d , m(0, •) = m 0 , φ(T, •) = φ T , in T d .
With a well-chosen time rescaling, one can introduce stationary MFG systems as long time average of time dependent ones. These take the form

(1.3)          λ + H(x, ∇φ) = f (x, m), in T d , -∇ • (D ξ H(x, ∇φ)m) = 0, in T d , T d m dx = 1, T d φ dx = 0, m ≥ 0, a.e. in T d ,
where an additional variable λ ∈ R appears in the first equation, which plays the role of the ergodic constant. We refer to [START_REF] Cardaliaguet | Long time average of mean field games with a nonlocal coupling[END_REF][START_REF] Cardaliaguet | Long time average of mean field games[END_REF] and [CG15, Section 4] for results on the limiting procedure and well-posedness of (1.3).

The case when the running cost f is monotone and regularizing (nonlocal) in the measure variable, is well understood in the literature. Already in the original works of J.-M. Lasry and P.-L. Lions (see also [START_REF] Lions | Cours au Collège de France[END_REF]) has been shown the well-posedness of both systems (1.2) and (1.3). By contrast, when f is a local function of m, the question of well-posedness of these systems is more subtle. Notice also that in the lack of sufficient regularity for the density variable m (when f is a local function of m), the control problem (1.1) is not meaningful, while the system (1.2) may still have solutions in a suitable weak sense.

In the case when f is a local function of m, PDE techniques, which are used successfully for second order systems (see for instance [START_REF] Gomes | Time-dependent mean-field games in the subquadratic case[END_REF][START_REF] Gomes | Time-dependent mean-field games in the superquadratic case[END_REF][START_REF] Gomes | Regularity theory for mean-field game systems[END_REF][START_REF] Porretta | Weak solutions to Fokker-Planck equations and mean field games[END_REF]), can no longer be used to show the well-posedness of the system (1.2). Also, in general one cannot hope for the existence of classical solutions of (1.2). Nevertheless, as it was already mentioned in [START_REF] Lasry | Jeux à champ moyen II. Horizon fini et contrôle optimal[END_REF] that when f is non-decreasing in its second variable and H convex in its second variable, at least formally, systems (1.2) and (1.3) can be seen as first order optimality conditions of two convex optimization problems in duality. In the first order case (or when a degenerate diffusion is present) these arguments were made rigorous in a series of papers (see [START_REF] Cardaliaguet | Weak solutions for first order mean field games with local coupling[END_REF][START_REF] Cardaliaguet | Mean field games systems of first order[END_REF][START_REF] Cardaliaguet | Second order mean field games with degenerate diffusion and local coupling[END_REF]), and the existence of weak solutions (understood in a suitable sense) was shown. An important and interesting question in this context is the regularity of solutions of (1.2) and (1.3).

A first result in this direction was obtained in [START_REF] Cardaliaguet | Sobolev regularity for the first order Hamilton-Jacobi equation[END_REF], which, it must be emphasized, provides first and foremost a regularity result for viscosity solutions of first-order Hamilton-Jacobi equations. Using techniques involving the inverse Hölder inequality, the authors obtained W 1,1 loc ((0, T ) × T d ) estimates and a.e. differentiability of φ, depending only on the summability of the right-hand side of the equation. In the context of mean field games, this amounts to requiring a sufficiently high growth condition on f (x, m). Notice that their estimates are not owing to the variational structure of the mean field games system; in particular the density variable plays no role.

A second direction in the search for Sobolev estimates for first order MFGs was initiated in [START_REF] Prosinski | Global-in-time regularity via duality for congestion-penalized mean field games[END_REF] (see also [START_REF] Santambrogio | Regularity via duality in calculus of variations and degenerate elliptic pdes[END_REF]), where the authors obtained H 1 loc ((0, T ] × T d ) estimates for a well-chosen function of m, where (φ, m) is a weak solution of (1.2). Their technique is based on a so-called regularity by duality method, which relies on the fact that (1.2) is the optimality condition for two convex optimization problems in duality. In fact, these techniques are inspired from [START_REF] Cardaliaguet | First order Mean Field Games with density constraints: pressure equals price[END_REF] and more precisely they have their roots in [START_REF] Brenier | Minimal geodesics on groups of volume-preserving maps and generalized solutions of the Euler equations[END_REF] (see also [START_REF] Ambrosio | On the regularity of the pressure field of Brenier's weak solutions to incompressible Euler equations[END_REF]), where Y. Brenier introduced a similar approach to obtain regularity estimates for the pressure field arising from generalized solutions of the incompressible Euler equations. In a nutshell, the basic idea of this method is the following. One considers special competitors for the optimization problem involving the m variable that are translations of the optimizers in time and space, then one compares the energy of these new competitors to the energy of the optimizers. Finally, assuming a sort of coercivity/monotonicity property of the running cost f in the second variable (see the precise assumption in Sections 2-3), one can obtain a differential quotient estimate for a function of m via the difference of the energies and deduce the Sobolev estimate. A drawback of this method is that since the optimization problem is subject to the continuity equation with an initial condition, when one constructs the new competitors (via the time-space translations), one must preserve the initial condition, so it seems impossible to construct translations that allow to vary the initial time as well. We shall give more details on this method in Section 2.

Let us mention that in the framework of first order models, P.-L. Lions in his lectures ( [START_REF] Lions | Cours au Collège de France[END_REF]) showed -at least formally -how to obtain a priori estimates in the case of the planning problem (where the initial and final density of the agents is prescribed). In particular, by rewriting the system as a very degenerate elliptic system in time and space, he obtains global L ∞ estimates for the density variable m and W 1,∞ estimates for the value function φ. However, his techniques require strong smoothness assumptions on the prescribed densities, Hamiltonian and coupling and positivity (and monotonicity) of the coupling function. Nevertheless, our objective in this paper is different: it is to obtain first order Sobolev estimates for the m variable and second order Sobolev estimates for φ. Thus, it seems that such an approach using degenerate elliptic equations to obtain these regularity estimates cannot be applied in our case. In the same context enters the recent paper [START_REF] Lavenant | Optimal density evolution with congestion: L ∞ bounds via flow interchange techniques and applications to variational mean field games[END_REF], where the authors obtain local in time L ∞ estimates for first order MFG models (with quadratic Hamiltonians), using variational techniques. Their approach, however, is very different from the one of Lions and from ours as well.

Actually, the same techniques as in [START_REF] Prosinski | Global-in-time regularity via duality for congestion-penalized mean field games[END_REF][START_REF] Santambrogio | Regularity via duality in calculus of variations and degenerate elliptic pdes[END_REF] can be used -more or less in a straight forward manner -for stationary problems (since there is no time involved there) as well, so as our preparatory results, we present how to gain Sobolev estimates for the solutions of (1.3). More precisely, we have Theorem 1.1 (Theorem 3.3, Theorem 3.4). We assume the standard regularity and growth conditions on the data f and H (see Sections 2 and 3 for the precise assumptions). Let (φ, m, λ) be a solution of (1.3). Then, one has a second order Sobolev estimate for φ, i.e.

J * (H(•, ∇φ) + λ) ∈ H 1 (T d ),
and a first order Sobolev estimate for m, i.e.

J(m) ∈ H 1 (T d ),
where J, J * : R → R are functions depending on f with a precise growth condition that depends on the growth condition of f in its second variable.

As our main result, we obtain global in time Sobolev estimates for weak solutions of first order MFG systems up to the initial (and final) time. While we keep the variational framework and rely on the convex duality, because of the obstruction described above, our approach is different from the one used in [START_REF] Prosinski | Global-in-time regularity via duality for congestion-penalized mean field games[END_REF]. Our main results can be summarized in an informal way in the theorem below.

Theorem 1.2 (Proposition 4.3, Corollary 4.5). We assume that m 0 , φ T ∈ W 2,∞ (T d ) and we assume some regularity and growth conditions on H (which has a growth like | • | r for some r > 1 in its second variable) and f (which has a growth like | • | q-1 for some q > 1 in its second variable). We refer to Sections 2 and 4 for the precise assumptions.

Then there exist a constant C > 0 depending only on the data such that

(i) m q 2 -1 ∇m L 2 ([0,T ]×T d ) ≤ C; (ii) m 1/2 D(j 1 (∇φ)) L 2 ([0,T ]×T d ) ≤ C; (iii) if r = 2, then ∂ t (m q/2 ) L 1 ([0,T ]×T d ) ≤ C.
Here j 1 : R d → R d is a function depending on H and has a growth like

| • | r/2 .
The core idea of our technique is the following. The space translations of the optimizers will actually solve some auxiliary optimization problems, where the data are the space translations of the data of the original problems. This observation leads us to obtain fully global in time differential quotient estimates in the space variable. Then using the continuity equation -when the Hamiltonian has a quadratic growth in its second variable -the first order Sobolev estimates in time for m, globally in time, will be a consequence.

Let us mention that while in [START_REF] Prosinski | Global-in-time regularity via duality for congestion-penalized mean field games[END_REF] the authors considered only the case when H(x, ξ) = 1 2 |ξ| 2 , most of the computations seem to be adaptable to more general Hamiltonians to obtain local in time Sobolev estimates, but to gain estimates up to the final time t = T , it seems crucial for them that D ξ H is a Lipschitz function. In comparison with this, our techniques allow us to consider quite general Hamiltonians (with no restriction on the growth condition) to obtain the global in time estimates for the space variable. Nevertheless, for the global in time estimates in the time variable, it seems that we have a similar obstruction.

Let us remark also that the estimate Theorem 1.2-(ii) seems to appear naturally in the context of Hamilton-Jacobi equations in connection with the nonlinear adjoint method introduced by L.C. Evans in [START_REF] Evans | Adjoint and compensated compactness methods for Hamilton-Jacobi PDE[END_REF]. Nevertheless, we believe that it is interesting to see this kind of second order estimates in the framework of first order Hamilton-Jacobi equations arising from MFGs, especially if we take into consideration that solutions to these equations a priori had to be understood in a very weak sense.

As a last remark, let us mention that it seems that our techniques to obtain global in time Sobolev estimates can be applied for other type of problems as well, in particular in the case of the first order planning problem (see [START_REF] Lions | Cours au Collège de France[END_REF]). These estimates -together with other results on this kind of systems -is the subject of a subsequent work ( [GM]).

The structure of the paper is simple. In Section 2 we collect the existing results on the wellposedness of first order variational MFG systems. The main assumptions, the convex duality results for the optimizations problems are also presented here. We also give a short summary on the main ideas of the technique "regularity by duality" presented in [START_REF] Prosinski | Global-in-time regularity via duality for congestion-penalized mean field games[END_REF] in the framework of MFG.

In Section 3 we present the Sobolev estimates for first order stationary MFG systems and we prove Theorem 1.1.

Finally, Section 4 is devoted to the proof of the global in time Sobolev estimates presented in Theorem 1.2.

Preliminary results and the regularity by duality method

We list the standing assumptions which are sufficient to ensure the well-posedness of both (1.2) and (1.3).

Standing assumptions: the following conditions on the data were used to get the result in [START_REF] Cardaliaguet | Mean field games systems of first order[END_REF].

(H1) (Conditions on the Hamiltonian) H :

T d × R d → R is continuous in both variables, convex
and differentiable in the second variable, with D ξ H continuous in both variables. Moreover, H has superlinear growth in the gradient variable: there exist r > 1 and C > 0 such that

(2.1) 1 rC |ξ| r -C ≤ H(x, ξ) ≤ C r |ξ| r + C, ∀ (x, ξ) ∈ T d × R d .
We denote by H * (x, •) the Fenchel conjugate of H(x, •), which, due to the above assumptions, satisfies

(2.2) 1 r ′ C |ζ| r ′ -C ≤ H * (x, ζ) ≤ C r ′ |ζ| r ′ + C, ∀ (x, ζ) ∈ T d × R d , where r ′ = r/(r -1) is the conjugate of r.
(H2) (Conditions on the coupling) Let f be continuous on T d × (0, ∞), strictly increasing in the second variable, and there exists a constant C > 0 satisfying

(2.3) 1 C |m| q-1 -C ≤ f (x, m) ≤ C|m| q-1 + C, ∀ m ≥ 1, ∀ x ∈ T d . (H3) (Conditions on the antiderivative of f ) We define F so that F (x, •) is an antiderivative of f (x, •) on (0, ∞), that is, (2.4) F (x, m) = m 1 f (x, s) ds, ∀ m > 0.
It follows that F is continuous on T d × (0, ∞), is strictly convex and differentiable in the second variable, and satisfies the growth condition

(2.5) 1 qC |m| q -C ≤ F (x, m) ≤ C q |m| q + C, ∀ m ≥ 1, ∀ x ∈ T d .
For m < 0 we set F (x, m) = +∞. We denote by F (x, 0) the limit lim m→0 + F (x, m), which may be finite or +∞. We will denote throughout the conjugate exponent of q by p = q ′ . Define F * (x, •) to be the Fenchel conjugate of F (x, •) for each x ∈ T d . Note that

(2.6) 1 pC |a| p -C ≤ F * (x, a) ≤ C p |a| p + C, ∀a ≥ 0, ∀ x ∈ T d .
2.1. Preliminary results on the stationary system. Let us recall here the stationary MFG system.

(2.7)

         (i) λ + H(x, Dφ) = f (x, m(x)), in T d , (ii) -∇ • (mD ξ H(x, ∇φ)) = 0, in T d , (iii)
T d m dx = 1; T d φ dx = 0; m ≥ 0, a.e in T d . Definition 2.1. We say that a triple (λ, φ, m) ∈ R × W 1,pr (T d ) × L q (T d ) is a weak solution of (2.7) if (i) m ≥ 0 a.e. in T d , T d m dx = 1, T d φ dx = 0 and mD ξ H(•, ∇φ) ∈ L 1 (T d ),
(ii) Equation (2.7)-(i) holds in the following sense:

(2.8) λ + H(x, ∇φ(x)) = f (x, m(x)) a.e. in {m > 0}

and

(2.9)

λ + H(x, ∇φ(x)) ≤ f (x, m) a.e. in T d , (iii) Equation (2.7)-(ii) holds: (2.10) -∇ • (mD ξ H(x, ∇φ)) = 0 in T d ,
in the sense of distributions.

In the sequel we summarize the results on the well-posedness of system (2.7). We define, on

K 0 := R × W 1,pr (T d ), the functional (2.11) A(λ, φ) = T d F * (x, λ + H(x, ∇φ(x))) dx -λ.
Our first optimization problem is

(2.12) inf (λ,φ)∈K 0 A(λ, φ)
To describe the second optimization problem, let us denote by K 1 the set of pairs (m, w)

∈ L 1 (T d )× L 1 (T d ; R d ) such that m ≥ 0 a.e. in T d , T d m(x) dx = 1
, and w is divergence free, i.e. it satisfies in the sense of distributions

(2.13) ∇ • w = 0, in T d .
We define on K 1 the functional

B(m, w) = T d m(x)H * x, - w(x) m(x) + F (x, m(x)) dx,
where we use the usual convention to define H * (x, -w/m), i.e.

(2.14) mH * x, -w m = +∞, if m = 0 and w = 0, 0, if m = 0 and w = 0.

Since H * and F are bounded below and m ≥ 0 a.e., the integral in B(m, w) is well defined in R ∪ {+∞}. The second optimal control problem is the following:

(2.15) inf (m,w)∈K 1 B(m, w) .

Lemma 2.2 ([CG15]

). We have

(2.16) min (λ,φ)∈K 0 A(λ, φ) = -min (m,w)∈K 1 B(m, w).
Moreover, the minimum in the right-hand side is achieved by a unique pair (m, w)

∈ K 1 satisfying (m, w) ∈ L q (T d ) × L r ′ q r ′ +q-1 (T d ).

Theorem 2.3 ([CG15]

). There exists at least one solution (λ, φ, m) to the stationary MFG system (2.7). Moreover, the pair (λ, m) is unique. If (m, w) ∈ K 1 is a minimizer of (2.15) and (λ, φ) ∈ K 0 is a minimizer of (2.12), then (λ, φ, m) is a solution of the MFG system (2.7) and w = -mD ξ H(•, ∇φ) a.e.. Conversely, any solution (λ, φ, m) of (2.7) is such that the pair (m, -mD ξ H(•, ∇φ)) is the minimizer of (2.15) while (λ, φ) is a minimizer of (2.11).

2.2.

Preliminary results on the time dependent system. Let us recall the time dependent MFG system that we will study.

(2.17)

   (i) -∂ t φ + H(x, ∇φ) = f (x, m), in (0, T ) × T d , (ii) ∂ t m -∇ • (mD ξ H(x, ∇φ)) = 0, in (0, T ) × T d , (iii) φ(T, x) = φ T (x), m(0, x) = m 0 (x), in T d .
First, we need to impose some additional standard assumption on the initial and final data.

(H4) (Conditions on the initial and final conditions) m 0 is a probability measure on T d which is absolutely continuous with respect to Lebesgue measure, having density (which we also call m 0 ) in C(T d ). We suppose moreover that m 0 > 0 on T d . We assume that φ T :

T d → R is a Lipschitz continuous function on T d . Definition 2.4. A pair (φ, m) ∈ BV ((0, T ) × T d ) × L q ((0, T ) × T d
) is called a weak solution to the system (2.17) if it satisfies the following conditions.

(1) ∇φ ∈ L r ([0, T ] × T d ) and the maps mf (•, m), mH * (•, -D ξ H(•, ∇φ)) and mD ξ H(•, ∇φ) are integrable, (2) φ satisfies a first-order Hamilton-Jacobi inequality

(2.18) -∂ t φ + H(x, ∇φ) ≤ f (x, m)
in the sense of distributions (tested against smooth non-negative functions), the boundary condition φ(T, •) ≤ φ T in the sense of trace and the following equality

(2.19) T 0 T d m(t, x) [H(x, ∇φ(t, x)) -∇φ(t, x) • D ξ H(x, ∇φ(t, x)) -f (x, m(t, x))] dx dt = T d (φ T (x)m(T, x)) -φ(0, x)m 0 (x)) dx
(3) m satisfies the continuity equation

(2.20) ∂ t m -∇ • (mD ξ H(x, ∇φ)) = 0 in (0, T ) × T d , m(0, •) = m 0 ,
in the sense of distributions.

We consider two optimal control problems which are in duality, see [START_REF] Cardaliaguet | Weak solutions for first order mean field games with local coupling[END_REF]. By the abuse of notion, we use the same notations A and B for the two functionals as in the stationary case. First, the control of the continuity equation reads as: minimize ). Second, we control the Hamilton-Jacobi equation: minimize

(2.21) B(m, w) = T 0 T d mH * x, - w m + F (x, m) dx dt + T d φ T m(T ) dx over K 1 = {(m, w) ∈ L 1 (T d ) × L 1 (T d ; R d ) : ∂ t m + ∇ • w = 0, m(0) = m 0 },
(2.22) A(φ) = T 0 T d F * (x, -∂ t φ + H(x, ∇φ)) dx dt - T d φ(0)m 0 dx over the set K 0 consisting of maps φ ∈ C 1 ([0, T ] × T d ) such that φ(T, •) = φ T . Set E 1 = C 1 ([0, T ] × T d ), E 0 = C([0, T ] × T d ; R) × C([0, T ] × T d ; R d ). On E 1 we define (2.23) F(φ) = - T d φ(0)m 0 + χ K 0 (φ). On E 0 define (2.24) G(a, b) = T 0 T d F * (x, -a(t, x) + H(x, b(t, x))) dx dt.
Then A can be written A(φ) = F(φ) + G(Λ(φ)) where Λ : E 1 → E 0 is given by Λ(φ) = (∂ t φ, ∇φ).

The two problems are in duality:

Theorem 2.5 (cf. [START_REF] Cardaliaguet | Weak solutions for first order mean field games with local coupling[END_REF]). We have

(2.25) inf φ∈K 0 A(φ) = -min (m,w)∈K 1 B(m, w).
Moreover, the minimum on the right-hand side is achieved by a unique pair (m, w)

∈ K 1 which must satisfy m ∈ L q ([0, T ] × T d ) and w ∈ L r ′ q r ′ +q-1 ([0, T ] × T d ).
Proof. Use the Fenchel-Rockafellar Theorem, which shows that

inf φ F(φ) + G(Λ(φ)) = max (m,w) -G * (-(m, w)) -F * (Λ * (m, w)).
The right-hand side is just -min B. See [START_REF] Cardaliaguet | Weak solutions for first order mean field games with local coupling[END_REF] for details.

We can relax the second problem in the following way. K will be defined as the set of all pairs (φ, α)

∈ BV ([0, T ] × T d ) × L 1 ([0, T ] × T d ) such that ∇φ ∈ L r ([0, T ] × T d ), φ(T, •) ≤ φ T in the sense of traces, α + ∈ L p ((0, T ) × T d ), φ ∈ L ∞ ([t, T ] × T d ) for every t ∈ (0, T ), and
(2.26)

-∂ t φ + H(x, ∇φ) ≤ α.
in the sense of distributions. Then on K we define the functional

A(φ, α) = T 0 T d F * (x, α(t, x)) dx dt - T d φ(0, x)m 0 (x) dx.
We have the following, due to [START_REF] Cardaliaguet | Mean field games systems of first order[END_REF].

Theorem 2.6. inf

φ∈K 0 A(φ) = min (φ,α)∈K
A(φ, α), the latter attained by some (φ, α) ∈ K of which φ can be obtained as a limit of smooth functions φ n .

Minimizers of A(φ, α) and B(m, w) precisely characterize weak solutions to the MFG system.

Theorem 2.7 (Existence and (partial) uniqueness of weak solutions). (i) If (m, w) ∈ K 1 is a minimizer of B and (φ, α) ∈ K is a minimizer of A, then (φ, m) is a weak solution of (2.17) and α(t, x) = f (x, m(t, x)) almost everywhere.

(ii) Conversely, if (φ, m) is a weak solution of (2.17), then there exist functions w, α such that (φ, α) ∈ K is a minimizer of A and (m, w) ∈ K 1 is a minimizer of B.

(iii) If (φ, m) and (φ ′ , m ′ ) are both weak solutions to (2.17), then m = m ′ almost everywhere while φ = φ ′ almost everywhere in the set {m > 0}.

Before we continue, let us remark that in [CG15], Theorem 2.7 holds only under the extra assumption r > d(q -1). The reason for this assumption was to have a priori upper bounds on the solution φ. However, in [START_REF] Cardaliaguet | Second order mean field games with degenerate diffusion and local coupling[END_REF] it was shown that this is not necessary; one has in general a priori bounds on the positive part φ + in L η for some η > 1. In particular, we have Theorem 2.8 (cf. Theorem 3.3 in [START_REF] Cardaliaguet | Second order mean field games with degenerate diffusion and local coupling[END_REF]). Suppose φ ∈ L r ((0, T ); W 1,r (T d )) satisfies, in the sense of distributions,

(2.27) (i) -∂ t φ + c 0 |∇φ| r ≤ α, in (0, T ) × T d , (ii) φ(T, x) ≤ φ T (x), a.e. in T d , for some r > 1, c 0 > 0, α ∈ L p ([0, T ] × T d ) and φ T ∈ L ∞ (T d ).
Then there exists a constant

C = C(p, d, r, c 0 , T, α L p ([0,T ]×T d ) , φ T L η (T d ) ) such that (2.28) φ + L ∞ ([0,T ];L η (T d )) + φ + L γ ([0,T ]×T d ) ≤ C
where η = d(r(p-1)+1) d-r(p-1) and γ = rp(1+d) d-r(p-1) if p < 1 + d r and η = γ = +∞ if p > 1 + d r . Using Theorem 2.8, one can modify the proof of Theorem 2.9 in [START_REF] Cardaliaguet | Mean field games systems of first order[END_REF] to get existence of minimizers for the relaxed problem min (φ,α)∈K A(φ, α). Indeed, let φ n ∈ K 0 be a minimizing sequence, so that (2.29)

T 0 T d F * (x, α n ) dx dt - T d φ n (0)m 0 dx → inf φ∈K 0 A(φ)
where

α n = -∂ t φ n + H(x, ∇φ n ). For ε ∈ R, set K ε = sup x∈T d F (x, ε
) and note that F * (x, a) ≥ aε -K ε for all a ∈ R. By the bounds (2.6) it follows that

(2.30)

T 0 T d F * (x, α n ) dx dt ≥ 1 pC (α n ) + p L p -ε T 0 T d (α n ) -dx dt -K ε T -C. Multiplying (α n ) -= (α n ) + + ∂ t φ n -H(x, ∇φ n
) by m 0 then integrating, using (2.1) and the fact that m 0 is bounded above and below by positive constants, we deduce

(2.31)

T 0 T d (α n ) -dx dt ≤ C T 0 T d (α n ) + m 0 dx dt + T d φ T m 0 -φ n (0)m 0 dx + CT ≤ C (α n ) + L p - T d φ n (0)m 0 dx + C.

Since

(2.32)

T 0 T d F * (x, α n ) dx dt - T d φ n (0)m 0 dx ≤ C
we can combine (2.30) and (2.31) to get

(2.33) 1 pC (α n ) + p L p -εC (α n ) + L p -(1 -ε) T d φ n (0)m 0 dx ≤ C + K ε T.
By fixing ε > 0 small enough and using (2.31) once more, we obtain a priori upper bounds on both (α n ) + p and -T d φ n (0)m 0 dx. On the other hand, since

(2.34) -

T d φ n (t)m 0 dx ≤ - T d φ n (0)m 0 dx + t 0 T d (α n -H(x, ∇φ n ))m 0 dx ≤ - T d φ n (0)m 0 dx + C (α n ) + p + C,
we also have an upper bound on -T d φ n (t)m 0 dx that holds uniformly in time. Combining this with Theorem 2.8, we deduce that

φ n is bounded in L ∞ ([0, T ]; L 1 (T d )). Moreover, (2.31) implies that (α n ) -is bounded in L 1 .
The rest of the proof now follows that of Theorem 2.9 in [START_REF] Cardaliaguet | Mean field games systems of first order[END_REF] and Proposition 5.4 in [START_REF] Cardaliaguet | Second order mean field games with degenerate diffusion and local coupling[END_REF]. In particular, a modified version of α n now has a weak limit α in L 1 ([0, T ] × T d ), and φ n is seen to have a weak limit φ in BV ([0, T ] × T d ), and one shows (cf. [START_REF] Cardaliaguet | Second order mean field games with degenerate diffusion and local coupling[END_REF]) that (φ, α) ∈ K is a minimizer of A.

The above arguments imply in particular that both optimization problems, and hence the MFG system have a solution, without imposing the joint assumption on the growth condition of f and H.

2.3.

The regularity by duality method. The idea to obtain Sobolev estimates by duality on the m variable in the system (1.2) used in [START_REF] Prosinski | Global-in-time regularity via duality for congestion-penalized mean field games[END_REF] dates back to [START_REF] Brenier | Minimal geodesics on groups of volume-preserving maps and generalized solutions of the Euler equations[END_REF] (see also [START_REF] Ambrosio | On the regularity of the pressure field of Brenier's weak solutions to incompressible Euler equations[END_REF] and [START_REF] Cardaliaguet | First order Mean Field Games with density constraints: pressure equals price[END_REF]) and it is as follows. Consider the two optimization problems involving the functionals A and B and suppose that these have minimizers (φ, α) and (m, w) respectively in the corresponding spaces as described above. The duality implies in particular that A(φ, α) + B(m, w) = 0. Now, consider η ∈ R small and δ ∈ T d such that |δ| is small and construct some special competitors by translations for the problem involving B, i.e. set (m η,δ , w η,δ )(t, x) := (m(t + η, x + δ), w(t + η, x + δ)).

Let us remark that in order to preserve in particular the initial condition m(0, •) = m 0 , one needs to use some special cutoff functions in time, when one considers time translations, but let us skip this detail in this informal description.

Suppose that F satisfies the strong coercivity condition (H6). Then, the key part of the analysis is to show that

(2.35) B(m η,δ , w η,δ ) ≤ B(m, w) + C(η 2 + |δ| 2 )
for some constant C > 0 independent of η and δ. Then, one can show some H 1 loc estimates for m. Indeed, first we will have

A( φ, α) + B( m, w) ≥ J( m) -J * (α) 2 L 2
for any competitors ( φ, α) and ( m, w) respectively. This means in particular by the duality that for the optimizers (φ, α) and (m, w) one has J(m) = J * (α). Second, by these last two observations if we add the quantity A(φ, α) to both sides of (2.35), one gets

J(m) -J(m η,δ ) 2 L 2 ≤ C(η 2 + |δ| 2
), thus the Sobolev estimate follows for J(m).

A warm-up: Sobolev estimates for the stationary problem

Our goal, in fact, is to show that the weak solutions (φ, m, λ) given by Theorem 2.3 are such that for φ one can derive a second order Sobolev estimate, while in the same time one can obtain first order Sobolev estimates for m. These results are in the spirit of those in [START_REF] Prosinski | Global-in-time regularity via duality for congestion-penalized mean field games[END_REF].

We will make some extra assumptions:

(H5) F * (•, α) ∈ C 1,1 (T d ), ∀ α ∈ R, D x F * (x, •) ∈ Lip loc (R), D α F * ∈ Lip loc (T d × R) and H(•, ξ) ∈ C 1,1 (T d ), ∀ ξ ∈ R d .
For h ∈ T d , we define the translates

φ h (x) := φ(x -h). Choose (λ, φ) ∈ K 0 such that A(λ, φ) is finite. Observe that A(λ, φ h ) = T d F * (x, λ + H(x, ∇φ(x -h))) dx -λ = T d F * (x + h, λ + H(x + h, ∇φ(x))) dx -λ is also finite. Lemma 3.1. Suppose (λ, φ) ∈ K 0 is such that A(λ, φ) is finite. Let A : T d → R be given by A(h) := A(λ, φ h ). Then A is C 1,1 (T d ).
Proof. Observe that

A(y + h) -A(y) = T d [F * (x + h, λ + H(x + h, ∇φ y (x))) -F * (x, λ + H(x, ∇φ y (x)))] dx = T d D x F * (x, λ + H(x, ∇φ y (x))) • h dx + T d D α F * (x, λ + H(x, ∇φ y (x)))D x H(x, ∇φ y (x)) • h dx + o(|h|). It follows that A ∈ C 1 (T d ) with DA(h) = T d [D x F * (x, λ + H(x, ∇φ h (x))) + D α F * (x, λ + H(x, ∇φ h (x)))D x H(x, ∇φ h (x))] dx = T d [D x F * (x+h, λ+H(x+h, ∇φ(x)))+D α F * (x+h, λ+H(x+h, ∇φ(x)))D x H(x+h, ∇φ(x))] dx.
By the assumption (H5) we see that DA(h) is Lipschitz continuous, as desired.

Lemma 3.2. Suppose (λ, φ) ∈ K 0 is such that A(λ, φ) is minimized. Then A(h) = A(λ, φ h ) = A(λ, φ) + O(|h| 2 ).
Proof. By the previous Lemma, we can expand A(h) = A(0) + DA(0)h + O(|h| 2 ), and since h = 0 is a minimizer for A we have DA(0) = 0. Now we are going to assume the following coercivity estimate on F : there exist J, J * : [0, +∞) → R and c 0 > 0 such that (H6)

F (x, m) + F * (x, a) -ma ≥ c 0 |J(m) -J * (a)| 2 .
Here the functions J, J * can be derived explicitly for many common examples of F (see the examples in [PS17, Section 3]). We deduce that

A(λ, φ) + B(m, w) = T d [F * (x, λ + H(x, ∇φ(x))) + m(x)H * (x, -w(x)/m(x)) + F (x, m(x))] dx -λ ≥ T d c 0 |J(m) -J * (λ + H(x, ∇φ(x)))| 2 dx + T d m(x) [H(x, ∇φ(x)) + H * (x, -w(x)/m(x))] dx ≥ T d c 0 |J(m) -J * (λ + H(x, ∇φ(x)))| 2 -∇φ(x) • w(x) dx = c 0 T d |J(m) -J * (λ + H(x, ∇φ(x)))| 2 dx for any (λ, φ) ∈ K 0 and (m, w) ∈ K 1 such that w ∈ L (q ′ r) ′ (T d ) = L r ′ q/(r ′ +q-1) (T d ) (note that ∇ • w = 0 in the sense of distributions, hence T d ∇φ(x) • w(x) dx = 0).
In particular, suppose that (λ, φ) ∈ K 0 and (m, w) ∈ K 1 are optimal, that is, A(λ, φ) and B(m, w) are minimized. Consider the translates φ h (x) := φ(x -h). Then we have the estimates

J * (H(•, ∇φ h ) + λ) -J * (H(•, ∇φ) + λ) 2 L 2 (T d ) ≤ 2 J * (H(•, ∇φ h ) + λ) -J(m) 2 L 2 (T d ) + 2 J(m) -J * (H(•, ∇φ) + λ) 2 L 2 (T d ) ≤ 2 c 0 (A(λ, φ h ) + B(m, w)) + 2 c 0 (A(λ, φ) + B(m, w)) = 2 c 0 (A(λ, φ h ) -A(λ, φ)) = O(|h| 2 ),
where in the last line we have used the duality result of Lemma 2.2 and the estimate from Lemma 3.2. This proves

Theorem 3.3. J * (H(•, ∇φ) + λ) is in H 1 (T d ).
Let us remark that the very same approach applies (even with simpler computations) to obtain Sobolev estimates for the m variable. Indeed, let us assume that the following regularity conditions hold true.

(H5') F (•, m) ∈ C 1,1 (T d ), ∀ m ∈ [0, +∞), and H * (•, ζ) ∈ C 1,1 (T d ), ∀ ζ ∈ R d .
Then, considering any competitor (m, w) such that B(m, w) is finite, for h ∈ T d defining the translates (m h , w h ), it is easy to show that the map

T d ∋ h → B(m h , w h ) is C 1,1 .
A similar reasoning as for the proof of Theorem 3.3 yields the estimate for J(m). More precisely, one has

Theorem 3.4. J(m) is in H 1 (T d ).
Observe also that formally Theorem 3.3 implies Theorem 3.4. In fact, since λ + H(x, ∇φ) = f (x, m), Theorem 3.3 implies that J * (f (•, m)) is in H 1 (T d ). Since the coupling f satisfying the growth condition (2.3), J grows like | • | q/2 and J * has a growth like | • | q ′ /2 , from where J * (f (•, m)) has exactly the same growth as J(m). However, this reasoning stays formal, one of the main obstructions being that λ + H(x, ∇φ) = f (x, m) is an equality a.e. only on {m > 0}. Thus, in order to obtain the Sobolev estimate for m, one needs to use the regularity by duality machinery, as described above.

Sobolev estimates for the time-dependent case

Our goal in this section is to prove global in time Sobolev estimates for the time dependent system (1.2). We will see that the technique is quite different from the duality method used in the previous section. The main reason for this is that we need to perturb the data, and not just the solution. In this context, let us assume some extra monotonicity/coercivity conditions for f and H.

Additional assumptions

(H7) (Conditions on the coupling) Let f be continuous on T d × (0, ∞), strictly increasing in the second variable, satisfying (2.3). Moreover, we will assume that f (x, m) is Lipschitz with respect to x, specifically

(4.1) |f (x, m) -f (y, m)| ≤ Cm q-1 |x -y| ∀x, y ∈ T d , m ≥ 0;
and also that f (x, m) is strongly monotone in m, i.e. there exists c 0 > 0 such that

(4.2) (f (x, m) -f (x, m)) ( m -m) ≥ c 0 min{ mq-2 , m q-2 }| m -m| 2 ∀ m, m ≥ 0, m = m.
(Note: if q < 2 one should interpret 0 q-2 as +∞ in (4.2). In this way, when m = 0, for instance, (4.2) reduces to f (x, m)m ≥ c 0 m q , as in the more regular case q ≥ 2.) (H8) (Coercivity assumptions.) There exist j 1 , j 2 : R d → R d such that

(4.3) H(x, ξ) + H * (x, ζ) -ξ • ζ ≥ c 0 |j 1 (ξ) -j 2 (ζ)| 2 .
In particular, and in light of our restriction (2.1), we will posit that j 1 (p) = |p| r/2-1 p and j 2 (p) = |p| r ′ /2-1 p.

Remark 4.1. The following "canonical" examples satisfy all of the hypotheses on the data:

f (x, m) = c 1 (x)m q-1 , H(x, ξ) = c 2 (x)|ξ| r ,
where c 1 , c 2 are continuous functions bounded below by a positive constant, and c 1 is Lipschitz. In order to satisfy the hypotheses of Proposition 4.3, we would also require c 2 to have W 2,∞ regularity.

Remark 4.2. Hypothesis (H7) is akin to (H6), in that the monotonicity estimates on f (x, •) are directly related to the convexity of F (x, •). In particular, if F (x, m) = 1 q c 1 (x)m q with c 1 as described in the previous remark, then (H7) is satisfied (because f (x, m) = c 1 (x)m q-1 ) while (H6) is satisfied with J(m) = m q/2 . Then the result of Proposition 4.3 below can be interpreted as giving an H 1 estimate on J(m).

4.1. Regularity in space. Let δ ∈ T d be given. Here we will consider perturbations of the data which are translations in the space variable. To wit, we use the superscript δ to denote translation by δ, i.e. α δ (t, x) = α(t, x + δ), F δ (x, m) = F (x + δ, m), and so on for all other data. Accordingly, we will write

A δ = A F δ ,m δ 0 , B δ = B H δ ,F δ ,φ δ T , K δ = K(φ δ T , H δ ), K δ 1 = K 1 (m δ 0 )
. Let (φ, α) ∈ K be a minimizer for A and (m, w) ∈ K 1 be a minimizer for B. To construct minimizers for A δ and B δ , we can use translations in space of (φ, α) and (m, w). Indeed, we have that (φ δ , α δ ) ∈ K δ and that A δ (φ δ , α δ ) = A(φ, α). One can then deduce that (φ δ , α δ ) is a minimizer for A δ , since any competitor can be translated backwards to get a competitor for (φ, α), which is a minimizer for A. In the same way (m δ , w δ )(t, x) = (m, w)(t, x + δ) is a minimizer for B δ .

We will use these facts in the following proof.

Proposition 4.3. Let m 0 , φ T ∈ W 2,∞ (T d ), and assume that H is twice continuously differentiable in x with

(H9) |D 2 x H(x, ξ)| ≤ C|ξ| r + C. Then m q 2 -1 ∇m L 2 ([0,T ]×T d ) ≤ C and m 1/2 D(j 1 (∇φ)) L 2 ([0,T ]×T d ) ≤ C. Proof.
Step 1. As above, let (m, w) ∈ K 1 (m 0 ) be a minimizer for B = B H,F,φ T . Initially we take a minimizing sequence of smooth functions

φ n ∈ K 0 (φ T ) such that A(φ n ) → inf A. Here A = A F,m 0 .
Denoting the translates as above, we see that φ δ n ∈ K 0 (φ δ T ) is also a minimizing sequence for

A δ = A F δ ,m δ 0 .
Following the argument in [START_REF] Cardaliaguet | Mean field games systems of first order[END_REF] we get that, up to a subsequence,

φ n → φ in L γ ([0, T ] × T d ) with 1 ≤ γ < d/(d -1), ∇φ n ⇀ ∇φ weakly in L r ([0, T ] × T d ; R d ) as n → +∞, and for α n := -∂ t φ n + H(x, ∇φ n ) we have αn := α n χ {αn≥-ln} → α in L 1 ([0, T ] × T d ) as n → +∞, where l n → ∞ as n → +∞; moreover, (φ, α) ∈ K is a minimizer of A.
In particular, for all M > 0 we can pass to a subsequence on which

α n χ {αn≥-M } → αχ {α≥-M } in L p ([0, T ] × T d ).
We now show that (up to a subsequence) (∇φ n ) converges weakly to ∇φ in L r m ([0, T ] × T d ; R d ). To see this, first use φ n as a test function in ∂ t m + ∇ • w = 0 to get (4.4)

T d φ T m(T ) -φ n (0)m 0 = T 0 T d (H(x, ∇φ n ) -α n )m + ∇φ n • w dx dt.
Using (2.1) we get

(4.5) 1 C T 0 T d |∇φ n | r m ≤ φ T ∞ + T d |φ n (0)|m 0 + T 0 T d (α n ) + m + |∇φ n |m w m dx dt,
then apply Young's inequality to get

(4.6) 1 C T 0 T d |∇φ n | r m ≤ φ T ∞ + T d |φ n (0)|m 0 + T 0 T d (α n ) + m + Cm w m r ′ dx dt
where, as usual, C takes on a larger value. As argued in [START_REF] Cardaliaguet | Second order mean field games with degenerate diffusion and local coupling[END_REF], lim sup n→∞ T d |φ n (0)|m 0 ≤ C.

Since (α n ) + is bounded in L p and m w m r ′
is integrable, this means that the right-hand side of (4.6) is bounded by a constant, as desired. Thus, we can assume that (up to a subsequence) (∇φ n ) converges weakly to some ξ in L r m ([0, T ] × T d ; R d ). We now claim that ξ = ∇φ, m-a.e. in [0, T ] × T d . That is, we claim that mξ = m∇φ a.e. It is enough to show that (4.7)

T 0 T d ψ • ∇φ n m → T 0 T d ψ • ∇φ m for all ψ ∈ C ∞ c ((0, T ) × T d ; R d ).
Let K > 0 (large) and ε > 0 (small). Since (∇φ n ) converges weakly to ∇φ in L r ([0, T ] × T d ; R d ), we have

(4.8) {m≤K} ψ • ∇φ n m dx dt → {m≤K} ψ • ∇φ m dx dt.
On the other hand, Young's inequality yields (4.9)

{m≥K} ψ • ∇φ n m dx dt ≤ ε r r ψ ∞ T 0 T d |∇φ n | r m dx dt + ε -r ′ r ′ ψ ∞ {m≥K}
m dx dt so that, by (4.6), we have (4.10) lim sup

n→∞ {m≥K} ψ • ∇φ n m dx dt ≤ C ε r r ψ ∞ + ε -r ′ r ′ ψ ∞ {m≥K} m dx dt. Since m is bounded in L q ([0, T ] × T d ), one has {m≥K} m dx dt ≤ {m≥K} 1 dx dt 1/q ′ m L q ,
thus Chebyshev's inequality yields that {m≥K} m dx dt → 0, as K → +∞. Thus, we let K → ∞ and then ε → 0 to get in (4.9) lim sup Step 2. Now use φ δ n and φ n as test functions in ∂ t m+∇•w = 0 and ∂ t m δ +∇•w δ = 0 respectively to get (4.11)

T d φ δ T m(T ) -φ δ n (0)m 0 = T 0 T d (H(x + δ, ∇φ δ n ) -α δ n )m + ∇φ δ n • w dx dt and
(4.12)

T d φ T m δ (T ) -φ n (0)m δ 0 = T 0 T d (H(x, ∇φ n ) -α n )m δ + ∇φ n • w δ dx
T d (φ δ T -φ T )m(T ) -(φ δ n (0) -φ(0))m 0 = T 0 T d (H(x + δ, ∇φ δ n ) + H * (x, -w/m) + ∇φ δ n • w/m -α δ n + f (m))m dx dt = T 0 T d (H(x, ∇φ δ n ) + H * (x, -w/m) + ∇φ δ n • w/m -α δ n + f (m))m dx dt + T 0 T d 1 0 D x H(x + sδ, ∇φ δ n ), δ ds dx dt.
Likewise, combine (4.12) with (2.19) (but applied to translates) to get (4.14)

T d (φ T -φ δ T )m δ (T ) -(φ n (0) -φ δ (0))m δ 0 = T 0 T d (H(x, ∇φ n ) + H * (x + δ, -w δ /m δ ) + ∇φ n • w δ /m δ -α n + f δ (m δ ))m δ dx dt = T 0 T d (H(x -δ, ∇φ -δ n ) + H * (x, -w/m) + ∇φ -δ n • w/m -α -δ n + f (m))m dx dt = T 0 T d (H(x, ∇φ -δ n ) + H * (x, -w/m) + ∇φ -δ n • w/m -α -δ n + f (m))m dx dt - T 0 T d 1 0 D x H(x -sδ, ∇φ -δ n ), δ ds dx dt.
Note that, by the changes of variables x → x + δ for the first integral and x → x -δ for the second, followed by the translation s → 1 -s, we get

T 0 T d 1 0 D x H(x + sδ, ∇φ δ n ), δ ds dx dt - T 0 T d 1 0 D x H(x -sδ, ∇φ -δ n ), δ ds dx dt = T 0 T d 1 0 D x H(x + sδ, ∇φ n ) -D x H(x -sδ, ∇φ n ), δ ds dx dt = T 0 T d 1 0 s -s D 2 xx H(x + rδ, ∇φ n )δ, δ dr ds dx dt.
Then adding together (4.13) and (4.14) we get (4.15)

T 0 T d H(x, ∇φ δ n ) + H(x, ∇φ -δ n ) + 2H * (x, -w/m) + ∇φ δ n • w/m + ∇φ -δ n • w/m m dx dt = T d φ δ T + φ -δ T -2φ T m(T ) dx - T d (φ n (0)(m δ 0 + m -δ 0 ) -2φ(0)m 0 ) + T 0 T d α δ n + α -δ n -2f (m) m dx dt + T 0 T d 1 0 s -s D 2 xx H(x + rδ, ∇φ n )δ,
T 0 T d H(x, ∇φ δ ) + H(x, ∇φ -δ ) m dx dt ≤ lim inf n→+∞ T 0 T d H(x, ∇φ δ n ) + H(x, ∇φ -δ n ) m dx dt.
Then, by the previous weak convergence and by the fact that w/m ∈ L r ′ m ([0, T ] × T d ; R d ), we have that (4.17)

T 0 T d ∇φ δ • w/m + ∇φ -δ • w/m m dx dt = lim n→+∞ T 0 T d ∇φ δ n • w/m + ∇φ -δ n • w/m m dx dt Second, let us compute (4.18) lim inf n→∞ T d φ δ T + φ -δ T -2φ T m(T ) dx - T d φ n (0)(m δ 0 + m -δ 0 ) -2φ(0)m 0 dx + T 0 T d α δ n + α -δ n -2f (m) m dx dt + T 0 T d 1 0 s -s D 2 xx H(x + rδ, ∇φ n )δ, δ dr ds dx dt . ≤ lim sup n→∞ T d φ δ T + φ -δ T -2φ T m(T ) dx - T d φ n (0)(m δ 0 + m -δ 0 ) -2φ(0)m 0 dx + T 0 T d α δ n + α -δ n -2f (m) m dx dt + T 0 T d 1 0 s -s D 2 xx H(x + rδ, ∇φ n )δ, δ dr ds dx dt . Let us recall that (φ n ) n≥0 converges in L 1 ([0, T ] × T d ) to φ ∈ BV ([0, T ] × T d ) and (φ n (0)) n≥0 is bounded in L 1 (T d ). In particular ∂ t φ n ⇀ ∂ t φ, as n → +∞, weakly-⋆ in M ([0, T ] × T d ). These imply furthermore that (4.19) lim sup n→+∞ - T d φ n (0)(m δ 0 + m -δ 0 ) dx ≤ - T d φ(0)(m δ 0 + m -δ 0 ) dx
Indeed, on the one hand we have As for the third term, we take an arbitrary M > 0 and get lim sup

- T d φ n (0, x)ψ(x) dx = - T d φ T (x)ψ(x) dx + T 0 T d ∂ t φ n ψ dx dt → - T d φ T (x)ψ(x) dx + T 0 T d ∂ t φψ dx dt = T d (φ(T, x) -φ T (x))ψ(x) dx - T d φ(0, x)ψ(x)
n→∞ T 0 T d α δ n + α -δ n -2f (m) m dx dt ≤ lim sup n→∞ T 0 T d α δ n χ {α δ n ≥-M } + α -δ n χ {α -δ n ≥-M } -2f (m) m dx dt = T 0 T d f δ (m δ )χ {f δ (m δ )≥-M } + f -δ (m -δ )χ {f -δ (m -δ )≥-M } -2f (m) m dx dt, which after letting M → ∞ becomes lim sup n→∞ T 0 T d α δ n + α -δ n -2f (m) m dx dt ≤ T 0 T d f δ (m δ ) + f -δ (m -δ ) -2f (m) m dx dt.
Lastly, by (H9) we can assert

T 0 T d |j 1 (∇φ δ ) -j 1 (∇φ -δ )| 2 m dx dt ≤ T d φ δ T + φ -δ T -2φ T m(T ) dx - T d φ(0) m δ 0 + m -δ 0 -2m 0 dx + T 0 T d f δ (m δ ) + f -δ (m -δ ) -2f (m) m dx dt + C|δ| 2 .
Step 3. We estimate the terms in (4.21). First, we have (4.22)

T d φ(0) 2m 0 -m δ 0 -m -δ 0 dx ≤ C|δ| 2 T d |φ(0)| dx,
where C depends on m 0 W 2,∞ , and we know that T d |φ(0)| dx is well-defined and finite (see [CGPT15, Lemma 5.1]). Similarly, (4.23)

T d m(T ) 2φ T -φ δ T -φ -δ T dx ≤ C|δ| 2 T d m(T ) dx = C|δ| 2 ,
where C depends on φ T W 2,∞ . There is one more term to estimate. We rewrite it using change of variables and then apply Equations (4.2) and (4.1) to get (4.24)

T 0 T d (f δ (m δ ) + f -δ (m -δ ) -2f (m))m dx dt = - T 0 T d (f δ (m δ ) -f (m))(m δ -m) dx dt = - {m δ ≤m} (f δ (m δ ) -f (m δ ))(m δ -m) dx dt - {m δ ≤m} (f (m δ ) -f (m))(m δ -m) dx dt - {m<m δ } (f δ (m δ ) -f δ (m))(m δ -m) dx dt - {m<m δ } (f δ (m) -f (m))(m δ -m) dx dt ≤ C T 0 T d |δ| min{(m δ ) q-1 , m q-1 }|m δ -m| dx dt -c 0 T 0 T d
min{(m δ ) q-2 , m q-2 }|m δ -m| 2 dx dt

≤ C|δ| 2 T 0 T d min{m δ , m} q dx dt - c 0 2 T 0 T d
min{(m δ ) q-2 , m q-2 }|m δ -m| 2 dx dt, where, we used Young's inequality in the last inequality, and the expression min{(m δ ) q-2 , m q-2 }|m δ -m| 2 is treated as zero whenever m δ = m (even in the case q < 2). Since min{(m δ ) q-2 , m q-2 }|m δ -m| 2 dx dt ≤ C|δ| 2 , where C depends only on the data. The result now follows from dividing both sides by |δ| 2 and letting δ → 0.

As a potential application of Proposition 4.3, one might hope to derive second order in space regularity for φ. For this, one would like some summability estimate on m -1 . Since, thanks to [START_REF] Cardaliaguet | Mean field games systems of first order[END_REF], the growth of f (x, m) is unrestricted for small values of m, all of our results are still valid for data satisfying, for example, F (x, m) ≥ 1 C m -s -C for some s ≥ 1, thus providing the desired estimate. Such a case would correspond to extreme congestion penalization, in which players experience a great benefit by moving toward unoccupied spaces. In this context, the following corollary is meaningful. m q/2-1 |∇ • w| ≤ Cm q/2-1 |∇m|(|∇φ| + 1) + Cm q/2 (|∇φ| + 1) + Cm q/2-1/2 m 1/2 |D 2 φ|.

Now since H(x, ∇φ) L 1 ≤ C, we have that |∇φ| ∈ L 2 ([0, T ] × T d ) by (2.1). Moreover, since q > 1 and m L q ≤ C it follows (by Hölder's inequality) that m q/2-1/2 L 2 ≤ C. Therefore (4.28) implies, by Hölder's inequality, that (4.31) m q/2-1 ∇ • w L 1 ([0,T ]×T d ) ≤ C.

We deduce the following: m q/2-1 ∇m L 2 , m 1/2 (|∇φ| r/2-1 D 2 φ + (r/2 -1)|∇φ| r/2-2 D 2 φ∇φ ⊗ ∇φ) L 2 ≤ C.

We compute (4.33) -∇ • w = ∇m • |∇φ| r-2 ∇φ + m|∇φ| r-2 ∆φ + m(r -2)|∇φ| r-4 (D 2 φ∇φ) • ∇φ.

If r < 2, then the last two terms of the r.h.s. are degenerate: it is possible that m|∇φ| r-2 ∆φ is not integrable, since even though one would have L 2 regularity for |∇φ| r/2-1 D 2 φ, there are no estimates on the remaining factor of m|∇φ| r/2-1 . On the other hand, if r > 2, then the first term on the r.h.s. is problematic: |∇φ| r-2 ∇φ ∈ L r ′ ([0, T ] × T d ), but r ′ < 2, and thus L 2 -regularity for ∇(m q/2 ) is not enough to ensure that the product ∇(m q/2 ) • |∇φ| r-2 ∇φ is integrable.

  where the continuity equation holds in the sense of distributions. (As usual, for (m, w) ∈ K 1 , the measure-valued function [0, T ] ∋ t → m(t) ∈ P(T d ) is continuous in time in the sense of weak convergence of measures, cf.[START_REF] Ambrosio | Existence, uniqueness, stability and differentiability properties of the flow associated to weakly differentiable vector fields[END_REF]

  n m dx dt = 0, In fact, by similar arguments, lim sup K→∞ lim sup n→∞ {m≥K} ψ • ∇φ n m dx dt -{m≥K} ψ • ∇φ m dx dt = 0, which, combined with (4.8), proves the claim.

  H(x + rδ, ∇φ n )δ, δ dr ds dx dt ≤ C( ∇φ n r r + 1)|δ| 2 ≤ C|δ| 2 . Now, the assumption (4.3) on H together with (4.16) and (4.17) on the one hand, the inequality |a + b| 2 ≤ 2(|a| 2 + |b| 2 ), ∀a, b ∈ R d and the previous estimations on the other hand yield (4.21) c 0 2

  δ , m} q dx dt ≤ T 0 T d m q dx dt ≤ C,we see from plugging (4.22), (4.23), and (4.24) into (4∇φ δ ) -j 1 (∇φ -δ )| 2 m dx dt

Corollary 4. 4 .(

 4 Suppose m -1 ∈ L s ([0, T ] × T d ) for some s ≥ 1. Then (4.26) ∇(j 1 (∇φ)) 2 2s s+1 ≤ C m -1 s ,where, for simplicity of writing • η denotes the norm• L η ([0,T ]×T d ) .Proof. Observe that by Hölder's inequality and Proposition 4.3, we have (4.27)j 1 (∇φ δ ) -j 1 (∇φ -δ ) |j 1 (∇φ δ ) -j 1 (∇φ -δ )| 2 )m dx dt ≤ C m -1 s |δ| 2 ,which implies the thesis of this corollary dividing both sides by |δ| 2 and taking δ → 0. 4.2. Some time regularity. Assume that r = 2, so that the Hamiltonian is of quadratic growth.In line with this assumption, we will assume that(H10) |D 2 pp H(x, ξ)| ≤ C, |D 2 xp H(x, ξ)| ≤ C(|ξ| + 1) ∀x ∈ T d , ∀ξ ∈ R d . Then the result of Proposition 4.3 reads (4.28) m q/2-1 ∇m L 2 ([0,T ]×T d ) , m 1/2 D 2 φ L 2 ([0,T ]×T d ) ≤ C.Recalling that w = -mD ξ H(x, ∇φ), we compute (4.29) -∇ • w = ∇m • D ξ H(x, ∇φ) + m tr(D 2 xp H(x, ∇φ)) + m tr(D 2 pp H(x, ∇φ)D 2 φ) and deduce from (H10) that (4.30)

  Corollary 4.5 (Global L 1 bounds on ∂ t m). Assume r = 2 and that (H10) holds. Suppose (m, w) is the minimizer of B. Then there exist C > 0 depending only on the data such that∂ t (m q/2 ) L 1 ([0,T ]×T d ) ≤ C.Proof. Apply (4.31) to the equation ∂ t m = -∇ • w. The argument breaks down if r = 2. For simplicity, set H(x, ξ) = 1 r |ξ| r . The result of Proposition 4.3 reads (4.32)

  δ dr ds dx dt. Now we pass to the limit as n → +∞ in each term above. First, we since H is convex in the second variable and (∇φ ±δ n ) converges weakly to ∇φ ±δ in L r m ([0, T ]×T d ; R d ), by weak lower semicontinuity we get (4.16)

  dx as n → +∞ for all ψ ∈ C(T d ). On the other hand, since φ(T, x) ≤ φ T (x) for a.e. x ∈ T d , for ψ ≥ 0 one can conclude that

	lim sup n→+∞	-	T d	φ n (0)ψ dx ≤ -	T d	φ(0)ψ dx,
	which implies (4.19) as desired.

Acknowledgements

We thank Filippo Santambrogio for the fruitful discussions during this project. In particular, some discussions with him and his paper [San17] inspired us in writing this paper. We thank also Pierre Cardaliaguet for his remarks and valuable suggestions on the manuscript.

The first author was supported by the National Science Foundation through Grant DMS-1612880. The second author was partially supported by the Gaspar Monge Program for Optimization and Operation Research (PGMO) via the project VarMFGPDE. He would like to thank the hospitality of Tohoku University and Tohoku Forum for Creativity, Sendai, Japan in the framework of the Thematic Program on Nonlinear Partial Differential Equations for Future Applications in Summer 2017.

Department of Mathematics, Baylor University, One Bear Place, Waco, TX 97328, USA E-mail address: jameson graber@baylor.edu Department of Mathematics, UCLA, 520 Portola Plaza, Los Angeles, CA 90095, USA E-mail address: alpar@math.ucla.edu