Albert C J Luo

Dennis O'connor

PERIODIC MOTIONS AND CHAOS WITH IMPACTING CHATTER AND STICK IN A GEAR TRANSMISSION SYSTEM

Keywords: Discontinuous dynamical systems, gear transmission, impacting chatter, stick motion, grazing bifurcation, periodic motions

] are adopted to determine the existence of periodic motions, which cannot be achieved from the local stability analysis. Numerical simulations are performed to illustrate periodic motions and stick motion criteria. Such an investigation may provide some clues to reduce the noise in gear transmission systems.

Introduction

In [START_REF] Luo | Nonlinear dynamics of a gear transmission dynamical system, Part I: Mechanism of impacting chatter and stick[END_REF], the motion mechanism of impacting chatter with stick in the gear transmission dynamical system was discussed. The reduction of noise in such a gear transmission system requires the identification of the existence of periodic motions. Hence, in this paper, periodic motions in such a gear transmission system will be addressed. The early investigation on periodic motions in gear transmission systems can be found in publications in the 1930's. den [START_REF] Den Hartog | Forced vibrations with non-linear spring constants[END_REF] gave a solution for symmetric, periodic motion of such a gear system. The used model is a piecewise system without damping. After half a century, [START_REF] Shaw | A periodically forced piecewise linear oscillator[END_REF] used the mapping technique to investigate periodic motion and stability of a piecewise linear system with a single discontinuity. [START_REF] Natsiavas | Periodic response and stability of oscillators with symmetric trilinear restoring force[END_REF] investigated the stability and bifurcation of periodic motions in a piecewise linear system with a symmetric trilinear spring. Based on the piecewise linear model, periodic motions for gear transmission systems were extensively investigated in 1990's (e.g. [START_REF] Comparin | Nonlinear frequency response characteristics of an impact pair[END_REF][START_REF] Kahraman | Nonlinear dynamics of a spur gear pair[END_REF][START_REF] Rook | Dynamic analysis of a reverse-idler gear pair with concurrent clearances[END_REF][START_REF] Theodossiades | Non-linear dynamics of gear-pair systems with periodic stiffness and backlash[END_REF]). The asymmetric motion of the piecewise linear system with the infinite stiffness was investigated numerically in [START_REF] Li | Bifurcation and Chaos in a forced zero-stiffness impact oscillator[END_REF]. The periodic solutions and bifurcations in the piecewise linear oscillator were further discussed, for which refer to [Klezka et al., 1992]. In recent years, the mapping concepts were extensively adopted to determine the periodic motion in impacting systems (e.g. [START_REF] Han | Chaotic motion of a horizontal impact pair[END_REF][START_REF] Luo | An unsymmetrical motion in a horizontal impact oscillator[END_REF]). This methodology was also employed to investigate the periodic and chaotic motions of the periodically driven piecewise linear systems (e.g. [START_REF] Luo | Global chaos in a periodically forced, linear system with a deadzone restoring force[END_REF][START_REF] Menon | A global period-1 motion of a periodically forced, piecewise linear system[END_REF]). In [START_REF] Luo | The mapping dynamics of periodic motions for a three-piecewise linear system under a periodic excitation[END_REF], a generalized methodology was given for predicting arbitrary periodic motions of discontinuous systems. A summary of the generalized method for the analytical prediction of complicated periodic motions in discontinuous systems was presented in [START_REF] Luo | Singularity and Dynamics on Discontinuous Vector Fields[END_REF]. [START_REF] Pfeiffer | Mechanische Systems mit unstetigen Ubergangen[END_REF] presented an impact model of gear transmissions, and the theoretical and experimental investigations on regular and chaotic motions in the gear box were carried out in [START_REF] Karagiannis | Theoretical and experimental investigations of gear box[END_REF]. [START_REF] Nordmark | Non-periodic motion caused by grazing incidence in an impact oscillator[END_REF] investigated the grazing-induced nonperiodic motion of the impacting oscillator, but only the necessary condition for such a grazing motion in the impacting oscillator was given. Furthermore, the normal form mapping for such grazing phenomena was further developed in [START_REF] Di Bernardo | Normal form maps for grazing bifurcation in n-dimensional piecewise-smooth dynamical systems[END_REF][START_REF] Di Bernardo | Bifurcation of dynamical systems with sliding: Derivation of normal form mappings[END_REF]. To model the gear transmission system, [START_REF] Luo | Periodic motion and grazing in a harmonically forced, piecewise, linear oscillator with impacts[END_REF] proposed a piecewise linear, impacting system to investigate the periodic motion. For such an oscillator, not only the simple periodic motion but also many complicated periodic motions exist. The switching between the two periodic motions causes motion complexity, and the asymmetric periodic motions exhibit the period-doubling bifurcation in such an oscillator. [START_REF] Luo | Arbitrary periodic motions and grazing switching of a forced piecewiselinear, impacting oscillator[END_REF] used the mapping technique to analytically predict arbitrary periodic motions for such an oscillator. As in [START_REF] Luo | Nonlinear dynamics of a gear transmission dynamical system, Part I: Mechanism of impacting chatter and stick[END_REF], it was stated that the impacting location in the model of [START_REF] Luo | Arbitrary periodic motions and grazing switching of a forced piecewiselinear, impacting oscillator[END_REF] was fixed and the perfect plastic impact was considered in this piecewise linear model. The separation of two gears occurs at the same location of gear impacts. To improve the modeling of gear transmission systems, the two gears will be considered to be independent, and impacts between the two separated oscillators will be adopted. Based on such a new model, the motion mechanism of the impacting chatter and stick was discussed in [START_REF] Luo | Nonlinear dynamics of a gear transmission dynamical system, Part I: Mechanism of impacting chatter and stick[END_REF]. For the further understanding of the origin of the vibration and noise in gear transmission systems, the corresponding periodic motions for such a new model will be investigated herein. Compared to the earliest results, the characteristics of complicated periodic motions can be predicted analytically by the mapping technique.

In this paper, basic mappings will be introduced first from the separation boundaries, and from the basic mappings, the mapping structures will be developed for any specific periodic motion. Further, the analytical prediction of periodic motions pertaining to impacting chatter and stick will be completed, and the corresponding local stability and bifurcation of the periodic motion will be analyzed. The grazing and stick conditions presented in [START_REF] Luo | Nonlinear dynamics of a gear transmission dynamical system, Part I: Mechanism of impacting chatter and stick[END_REF] will be employed for the existence of periodic motions. Numerical simulations will be carried out for illustration of periodic motions and stick criteria.

Basic Mappings

For the gear transmission system, equations of motion in the absolute frame are from [START_REF] Luo | Nonlinear dynamics of a gear transmission dynamical system, Part I: Mechanism of impacting chatter and stick[END_REF]

ẋ(i) α = F (i) α (x (i) α , t) ( 1 )
for i = 1, 2 and α = 1, 2, 3 with the following vectors

x (i) α = (x (i) α , ẋ(i) α ) T = (x (i) α , y (i) α ) T , F (i) α = ( ẋ(i) α , F (i) α ) T = (y (i) α , F (i) α ) T ;
(2)

where

F (i) α = -2ζ (i) α ẋ(i) α -(ω (i) α ) 2 x (i) α + b (i) α + Q (i) α cos Ωt, (3)
and the superscript "i" represents the ith mass and the subscript "α" represents the α-domain. Similarly, the corresponding equations of motion in the relative frame for i = 1, 2 and α = 1, 2, 3 are

ż(i) α = g (i) α (z (i) α , x ( ī) α , t), ẋ( ī) α = F ( ī) α (x ( ī) α , t) (4)
where i, ī ∈ {1, 2} and i = ī, and the relative vectors are

z (i) α = (z (i) α , ż(i) α ) T = (z (i) α , v (i) α ) T , g (i) α = ( ż(i) α , g (i) α ) T = (v (i) α , g (i) α ) T (5)
where

g (i) α = -2ζ (i) α ż(i) α -(ω (i) α ) 2 z (i) α + b (i) α + Q (i) α cos Ωt -ẍ( ī) α -2ζ (i) α ẋ( ī) α -(ω (i) α ) 2 x ( ī) α . ( 6
)
To investigate periodic impacting chatter with/ without stick in the gear transmission system, the mapping structure will be developed from separation boundaries in [START_REF] Luo | Nonlinear dynamics of a gear transmission dynamical system, Part I: Mechanism of impacting chatter and stick[END_REF]. Before the mapping structure for a prescribed impacting chatter motion is developed, the switching planes will be defined first, and from the switching planes, the basic mappings can be developed for the gear transmission system. From discontinuous boundaries in [START_REF] Luo | Nonlinear dynamics of a gear transmission dynamical system, Part I: Mechanism of impacting chatter and stick[END_REF], the switching planes based on the two impacting chatter boundaries are defined as

R Σ (i) 2∞ =      (t k , x (i) k , ẋ(i) k , ẋ( ī) k ) x ( ī) k = x (i) k - d 2 ẋ(i) k = ẋ( ī) k      , L Σ (i) 2∞ =      (t k , x (i) k , ẋ(i) k , ẋ( ī) k ) x ( ī) k = x (i) k + d 2 ẋ(i) k = ẋ( ī) k      .                      (7)
From now on, x

(i) k ≡ x (i) (t k ) and ẋ(i) k ≡ ẋ(i) (t k
) on the separation boundary at time t k are switching displacement and velocity in the absolute frame. The switching phase is defined by ϕ k = mod(Ωt k , 2π). In the relative frame, the switching planes are expressed as

R Σ (i) 2∞ = (t k , ż(i) k , x ( ī) k , ẋ( ī) k ) z (i) k = d 2 , ż(i) k = 0 L Σ (i) 2∞ = (t k , ż(i) k , x ( ī) k , ẋ( ī) k ) z (i) k = - d 2 , ż(i) k = 0          (8)
The two switching sets are then decomposed as

R Σ (i) 2∞ = R + Σ (i) 2∞ ∪ R -Σ (i) 2∞ and L Σ (i) 2∞ = L + Σ (i) 2∞ ∪ L -Σ (i) 2∞ (9)
where

R + Σ (i) 2∞ = (t k , ż(i) k , x ( ī) k , ẋ( ī) k ) z (i) k = d 2 , ż(i) k > 0 , R -Σ (i) 2∞ = (t k , ż(i) k , x ( ī) k , ẋ( ī) k ) z (i) k = d 2 , ż(i) k < 0 ; R + Σ (i) 2∞ = (t k , ż(i) k , x ( ī) k , ẋ( ī) k ) z (i) k = - d 2 , ż(i) k > 0 , R -Σ (i) 2∞ = (t k , z (i) k , x ( ī) k , ẋ( ī) k ) z (i) k = - d 2 , ż(i) k < 0 .                              (10)
Based on the above definitions of switching planes, four mappings are defined in the absolute frame as

P 2 : R Σ (i) 2∞ → R Σ (i) 2∞ , P 3 : R Σ (i) 2∞ → L Σ (i) 2∞ ; P 5 : L Σ (i) 2∞ → L Σ (i) 2∞ , P 6 : L Σ (i) 2∞ → R Σ (i) 2∞    (11)
and in the relative frame, the mappings are also defined as

P 2 : R -Σ (i) 2∞ → R + Σ (i) 2∞ , P 3 : R -Σ (i) 2∞ → L -Σ (i)
2∞

;

P 5 : L + Σ (i) 2∞ → L -Σ (i) 2∞ , P 6 : L + Σ (i) 2∞ → R + Σ (i) 2∞ .    (12)
Among four basic mappings, the two mappings (P 2 and P 4 ) are local and the other two mappings (P 3 and P 6 ) are global. The local mapping will map the motion from a switching plane onto itself. However, the global mapping will map the motion from a switching plane to another one. Such mappings are sketched in Fig. 1. The corresponding switching planes are labeled. The mappings for the absolute and relative frames are arranged in Figs. 1(a) and 1(b), respectively. On the impacting chatter boundaries, impacts are expressed by thin straight lines with arrows.

( ) i x ( ) i x ( ) 2 R i ∞ Σ ( ) 2 R i ∞ Σ ( ) 2 R i ∞ Σ ( ) 2 L i ∞ Σ ( ) 2 L i ∞ Σ ( ) 2 L i ∞ Σ 3 P 5 P 6 P 2 P ( ) i z ( ) i z ( ) 2 L i - ∞ Σ ( ) 2 R i - ∞ Σ ( ) 2 R i + ∞ Σ ( ) 2 L i + ∞
To investigate stick motions in the gear transmission system, the switching planes for stick are defined as

Σ (i) 12 = (t k , x (i) k , ẋ(i) k , ẋ( ī) k ) x (i) k = R x ( ī) 2-+ d 2 , ẋ(i) k = R ẋ( ī) 2-, Σ (i) 21 = (t k , x (i) k , ẋ(i) k , ẋ( ī) k ) x (i) k = R x ( ī) 2+ + d 2 , ẋ(i) k = R ẋ( ī) 2+ ; Σ (i) 23 = (t k , x (i) k , ẋ(i) k , ẋ( ī) k ) x (i) k = L x ( ī) 2-+ d 2 , ẋ(i) k = L ẋ( ī) 2-, Σ (i) 32 = (t k , x (i) k , ẋ(i) k , ẋ( ī) k ) x (i) k = L x ( ī) 2+ + d 2 , ẋ(i) k = L ẋ( ī) 2+                              (13)
for the absolute frame and

Σ (i) 12 = (t k , ż(i) k , x ( ī) k , ẋ( ī) k ) z (i) k = d 2 , ż(i) k = 0 -, Σ (i) 21 = (t k , ż(i) k , x ( ī) k , ẋ( ī) k ) z (i) k = d 2 , ż(i) k = 0 + ; Σ (i) 23 = (t k , ż(i) k , x ( ī) k , ẋ( ī) k ) z (i) k = - d 2 , ż(i) k = 0 -, Σ (i) 32 = (t k , ż(i) k , x ( ī) k , ẋ( ī) k ) z (i) k = - d 2 , ż(i) k = 0 +                              (14)
for the relative frame.

In the relative frame, the switching planes in Eqs. ( 10) and ( 13) are almost same. In Eq. ( 10), the relative velocity is nonzero (i.e. ż(i) k = 0) but in Eq. ( 13), the relative velocity is zero (i.e. ż(i) k = 0). The two switching planes can be treated as the same for all mappings. Except for two stick mappings (i.e. P 1 and P 4 ), the other mappings are the same as in Eqs. ( 11) and (12). From the stick switching planes, the mappings are defined as

P 1 : Σ (i) 21 → Σ (i) 12 , P 2 : Σ (i) 12 → Σ (i) 21 , P 3 : Σ (i) 12 → Σ (i) 23 ; P 4 : Σ (i) 23 → Σ (i) 32 , P 5 : Σ (i) 23 → Σ (i) 32 , P 6 : Σ (i) 32 → Σ (i) 21 .    ( 15
)
With mixed switching planes, six mappings are defined by

P 2 : Σ (i) 12 → R Σ (i) 2∞ , P 2 : R Σ (i) 2∞ → Σ (i) 21 ; P 3 : Σ (i) 12 → L Σ (i) 2∞ , P 3 : R Σ (i) 2∞ → Σ (i) 23 ; P 5 : Σ (i) 23 → L Σ (i) 2∞ , P 5 : L Σ (i) 2∞ → Σ (i) 32 , P 6 : Σ (i) 32 → R Σ (i) 2∞ , P 6 : L Σ (i) 2∞ → Σ (i) 21 .                  (16)
Similarly, the mappings based on the relative switching planes can be defined Set the vectors as ( ) i

y k ≡ (t k , x (i) k , ẋ(i) k , ẋ( ī) k ) T , w k ≡ (t k , ż(i) k , x ( ī) k , ẋ( ī) k ) T .    (17) ( ) i x ( ) i x ( ) 21 i Σ ( ) 12 i Σ ( ) 32 i Σ ( ) 23 i Σ 3 P 6 P 4 P 1 P 2 P 5 P ( ) i z ( ) i z ( ) L i Σ ( ) R i Σ ( ) 2 R i + ∞ 2 - ∞ 2 - ∞ Σ ( ) 2 L i + ∞ Σ 6 P 3 P 2 P 4 P ( ) 32 i Σ ( ) 23 i Σ ( ) 21 i Σ ( ) 12 i Σ (a) (b)
x For the impacting maps P σ (σ = 2, 3, 5, 6) in the absolute coordinate, y k+1 = P σ y k can be expressed by

( ) i x ( ) 21 i Σ ( ) 2 R i ∞ Σ ( ) 12 i Σ ( ) 32 i Σ ( ) 2 L i ∞ Σ ( ) 23 i Σ 3 P 5 P 6 P 2 P 2 P 5 P ( ) i x ( ) i x ( ) 21 i Σ ( ) 2 L i ∞ Σ ( ) 12 i Σ ( ) 32 i Σ ( ) 2 R i Σ ( ) 23
P σ : (t k , x (i) k , ẋ(i) k , ẋ( ī) k ) → (t k+1 , x (i) k+1 , ẋ(i) k+1 , ẋ( ī) k+1 ) (18)
From Appendix, the absolute displacement and velocity for two gear oscillators can be obtained with initial conditions (t k , x

(i) k , ẋ(i) k ) and (t k , x ( ī) k , ẋ( ī) k ).
The final state for time t k+1 can be given. The switching planes give x

( ī) γ = x (i) γ ± (d/2) (γ = k, k + 1)
. The four equations of displacement and velocity for two gears give a set of four algebraic equations, i.e.

f (σ) (y k , y k+1 ) = 0 (19)
where

f (σ) = (f (σ) 1 , f (σ) 2 , f (σ) 3 , f (σ) 4 ) T . ( 20
)
For the stick motion, the displacement and velocity of the ith-gear will be adopted. In addition, the following equations are used.

ẋ( ī) k+1 = ẋ(i) k+1 , g (i) α (x ( ī) k+1 , ẋ( ī) k+1 , ẍ( ī) k+1 , t k+1 ) = 0. (21)
With the condition x

( ī) γ = x (i) γ ±(d/2) (γ = k, k +1
), the algebraic equations in Eq. ( 19) can be obtained. If a mapping only starts or ends at the stick boundary, the corresponding displacement plus the following equation can be employed to obtain Eq. ( 19).

ẋ( ī) γ = ẋ(i) γ , for γ ∈ {k, k + 1}. ( 22
)
Based on the relative coordinate (z (i) , ż(i) ), the relative displacement and velocity can be used with the initial condition (t k , z

(i) k , ż(i) k )
on the switching boundary. The displacement and velocity with an initial condition (t k , x

( ī) k , ẋ( ī) k ) can be given. With conditions z (i) γ = ± (d/2) (γ = k, k + 1)
, the relative and absolute displacements and velocities generate a set of four algebraic equations as

h (σ) (w k , w k+1 ) = 0 (23)
where

h (σ) = (h (σ) 1 , h (σ) 2 , h (σ) 3 , h (σ) 4 ) T . ( 24
)
In a similar fashion, for stick mapping, we have

ż(i) k+1 = ż(i) k = 0, g (i) α (x ( ī) k+1 , ẋ( ī) k+1 , ẍ( ī) k+1 , t k+1 ) = 0.    (25)
If a mapping only starts or ends at the stick boundary, the relative displacement plus the following equation can be used to obtain Eq. ( 23).

ż(i) γ = 0, for γ ∈ {k, k + 1}. (26)
Finally, from Luo and O'Connor (1), the impact mapping on the impact boundaries is defined as

P 0 : R Σ (i) 2∞ → R Σ (i) 2∞ and P 0 : L Σ (i) 2∞ → L Σ (i) 2∞ (27)
in the absolute frame and

P 0 : R + Σ (i) 2∞ → R -Σ (i) 2∞ and P 0 : L -Σ (i) 2∞ → L + Σ (i) 2∞ (28)
in the relative frame. The corresponding functions in Eqs. ( 19) and ( 23), respectively, are

f (0) 1 = t k+1 -t k , f (0) 2 = x (i) k+1 -x (i) k , f (0) 3 = ẋ(i) k+1 -I (i) 1 ẋ(i) k -I (i) 2 ẋ( ī) k , f (0) 4 = ẋ( ī) k+1 -I ( ī) 1 ẋ( ī) k -I ( ī) 2 ẋ(i) k .                (29) h (0) 1 = t k+1 -t k , h (0) 2 = ż(i) k+1 + e ż(i) k , h (0) 3 = x ( ī) k+1 -x ( ī) k , h (0) 4 = ẋ( ī) k+1 - ẋ( ī) k - m i m 1 + m 2 ( ż(i) k+1 - ż(i) k ).                  (30)
For simplicity of mapping structures of periodic motions, the impact mapping will be dropped from now on, but the impact relation will be embedded.

Mapping Structures

To investigate periodic motions in such a gear transmission system, the notation for mapping actions of basic mappings is introduced as in [START_REF] Luo | The mapping dynamics of periodic motions for a three-piecewise linear system under a periodic excitation[END_REF][START_REF] Luo | Singularity and Dynamics on Discontinuous Vector Fields[END_REF]]

P n k •••n 2 n 1 ≡ P n k • • • • • P n 2 • P n 1 (31)
where the mapping P n j (n j ∈ {1, 2, . . . , 6}, j = 1, 2, . . . , k) is defined in the previous section. Consider a generalized mapping structure as

P (65 k s4 4 k s3 31 k s2 2 k s1 )•••(65 k 14 4 k 13 31 k 12 2 k 11 ) s-terms = P (65 k s4 4 k s3 31 k s2 2 k s1 ) • • • • • P (65 k 14 4 k 13 31 k 12 2 k 11 ) s-terms (32)
where (k µν ∈ {0, N}, µ = 1, 2, . . . , s, ν = 1, 2, 3, 4).

From the generalized mapping structure, consider a simple mapping structure of periodic motions for impacting chatter. For instance, the mapping structure is

P 65 n 32 m = P 6 • P 5 n • P 3 • P 2 m (33)
where m, n ∈ {0, N}. Such a mapping structure gives (m + 1)-impacts in boundary R ∂Ω (i) 2∞

and (n + 1)-impacts in boundary L ∂Ω (i) 2∞ , which are described by mappings P 2 and P 5 , respectively. Through the global mappings P 3 and P 6 , the impacting chatters on the two boundaries are connected together. Consider a periodic motion of P 65 n 32 m with period T 1 = k 1 T (k 1 ∈ N). If the mapping structure copies itself, the new mapping structure is:

P (65 n 32 m ) 2 = P 65 n 32 m • P 65 n 32 m (34)
The periodic motion of P (65 n 32 m ) 2 is obtained during a period of 2T 1 . In an alike fashion, such an action of mapping structure continues to copy itself with period-2 l T 1 .

P (65 n 32 m ) 2 l = P (65 n 32 m ) 2 l-1 • P (65 n 32 m ) 2 l-1 (35)
As l → ∞, a chaotic motion relative to mapping structure P 65 n 32 m is formed. The prescribed chaos is generated by period-doubling. However, if the grazing bifurcation occurs, such a mapping structure may not be copied by itself. The new mapping structures are combined by the two different mapping structures. For instance,

P 65 n 2 32 m 2 65 n 1 32 m 1 = P 65 n 2 32 m 2 • P 65 n 1 32 m 1 , . . . P 65 n l 32 m l •••65 n 1 32 m 1 = P 65 n l 32 m l • • • • • P 65 n 1 32 m 1 l-terms . ( 36
)
Such a gazing bifurcation will cause the discontinuity of periodic motions, and chaotic motions may exist between periodic motions of

P 65 n l 32 m l •••65 n 1 32 m 1 and P 65 n l-1 32 m l-1 •••65 n 1 32 m 1 .
For low excitation frequency, the impacting chatter accompanying stick motion exists in the gear transmission system. Consider a simple chatter with stick motion with the following mapping structure

P 645 n 312 m = P 6 • P 4 • P 5 n • P 3 • P 1 • P 2 m (37)
From the above mapping structure, m-impacts in the boundary R ∂Ω (i) 2∞ and n-impacts in the boundary L ∂Ω (i) 2∞ , are described by mappings P 2 and P 5 , respectively. In addition, both the mth mapping of P 2 and the nth mapping of P 5 map the impacting boundary to the stick boundary, and the corresponding stick mappings are P 1 and P 4 , respectively. The two global mappings P 3 and P 6 connect the impact and stick boundaries. Similarly, for period-doubling, the mapping structures are given by P (645 n 312 m ) 2 = P 645 n 312 m • P 645 n 312 m , . . .

P (645 n 312 m ) 2 l = P (645 n 312 m ) 2 l-1 • P (645 n 312 m ) 2 l-1 .      (38)
Due to grazing bifurcation, the mapping structures are:

P 645 n 2 32 m 2 65 n 1 312 m 1 = P 645 n 2 312 m 2 • P 645 n 1 312 m 1 , . . . P 645 n l 312 m l •••645 n 1 312 m 1 = P 645 n l 312 m l • • • • • P 645 n 1 312 m 1 l-terms . ( 39
)
To help one understand two sorts of mapping structures, the two simple mapping structures are shown in Figs. 4(a) and 4(b) for the impacting chatter with and without stick of two gear systems. Similarly, the other mapping structures can be discussed through the generalized mapping structure in Eq. (32). Further, the periodic and chaotic motions relative to a certain mapping structure can be determined.

( ) i x ( ) i x ••• 6 P ( ) 2 R i ∞ Σ ( ) 2 R i ∞ Σ ( ) 2 R i ∞ Σ ( ) 2 L i ∞ Σ ( ) 2 L i ∞ Σ ( ) 2 L i ∞ Σ 2 P 2 P 3 P 5 P 5 P ••• ••• ••• ( ) i x ( ) i x ••• ( ) 2 R i ∞ Σ ( ) 2 R i ∞ Σ ( ) 12 i Σ ( ) 32 i Σ ( ) 2 L i ∞ Σ ( ) 2 L i ∞ Σ ( ) 21 i Σ ( ) 23

Bifurcation Scenario

Before discussing impacting chatter in the gear transmission system, a bifurcation scenario is presented through the switching displacements, velocities and phases of two gear oscillators. All the numerical computations are completed from the closed-form solutions in Appendix. The parameters (m 1 = 2, m 2 = 1, r 1 = r 2 = 0.6, k 1 = 30, k 2 = 20, Q 0 = 30, d = 1.0 and e = 0.7) are considered. The impacting chatter with and without stick varying with excitation frequency are presented in Figs. 5(a)-5(d). The "C-chatter" and "S-chatter" represent the complicated and simple impacting chatters, respectively. The "Chatter w/stick" denotes the impacting chatter with stick for periodic motions. For high excitation frequency, it is observed that the impacting chatter motion possesses simple mapping structure, and the motion is very simple. The excitation frequency lies in the interval of Ω ∈ (4.159, 16.433) for the simplest, impacting-chatter, periodic motion relative to P 63 . However, the complex, impacting-chatter motion exists for Ω ∈ (1.398, 4.159). For Ω ∈ (0, 1.398), the impacting chatter with stick exists. The details of mapping structures for impacting chatter with stick are tabulated in Table 1, and the excitation frequency range for impacting chatter is in Table 2.

Periodic motion and stability

From mapping structures of periodic motions, the switching sets for any periodic motion can be determined through solving a set of nonlinear algebraic equations. Consider a periodic motion of mapping structure P (65 k s4 4 k s3 31 k s2 2 k s1 )•••(65 k 14 4 k 13 31 k 12 2 k 11 ) and the following relation holds.

P (65 k s4 4 k s3 31 k s2 2 k s1 )•••(65 k 14 4 k 13 31 k 12 2 k 11 ) y k = y k+2s+ P s m=1 P 4 j=1 k mj ( 40
)
where

y k = (t k , x (i) k , ẋ(i) k , ẋ( ī) k ) T . A set of vector equations is as f (2) (y k+1 , y k ) = 0, f (2) (y k+2 , y k+1 ) = 0, . . . f (6) (y k+2s+ P s m=1 P 4 j=1 k mj , y k+2s+ P s m=1 P 4 j=1 k mj -1 ) = 0;                    (41)
where

f (σ) = (f (σ) 1 , f (σ) 2 , f (σ) 3 , f (σ)
4 ) T is relative to governing equations of mapping P σ (σ ∈ {1, 2, . . . , 6}). The periodicity of the period-1 motion per N -periods requires

y k+2s+ P s m=1 P 4 j=1 k mj = y k (42) or x (i) k+2s+ P s m=1 P 4 j=1 k mj = x (i) k , ẋ(i) k+2s+ P s m=1 P 4 j=1 k mj = ẋ(i) k , ẋ( ī) k+2s+ P s m=1 P 4 j=1 k mj = ẋ( ī) k , Ωt k+2s+ P s m=1 P 4 j=1 k mj ≡ Ωt k + 2Nπ.                  (43)
Solving Eqs. ( 41) and ( 42) generates the switching sets for periodic motions. Once the analytical prediction of any periodic motion is obtained, the corresponding stability and bifurcation analysis can be completed.

The local stability and bifurcation for such a period-1 motion is determined through the corresponding Jacobian matrix of the Poincaré mapping. From Eq. ( 40), the Jacobian matrix is computed by the chain rule, i.e.

DP = DP

(65 k n4 4 k n3 31 k n2 2 k n1 )•••(65 k 14 4 k 13 31 k 12 2 k 11 ) = n m=1 (DP 6 • DP k m4 5 • DP k m3 4 • DP 3 • DP k m2 1 • DP k m1 1 ) (44)
where

DP λ =                   ∂t ν+1 ∂t v ∂t ν+1 ∂x (i) ν ∂t ν+1 ∂y (i) ν ∂t ν+1 ∂y ( ī) ν ∂x (i) ν+1 ∂t v ∂x (i) ν+1 ∂x (i) ν ∂x (i) ν+1 ∂y (i) ν ∂x (i) ν+1 ∂y ( ī) ν ∂y (i) ν+1 ∂t v ∂y (i) ν+1 ∂x (i) ν ∂y (i) ν+1 ∂y (i) ν ∂y (i) ν+1 ∂y ( ī) ν ∂y ( ī) ν+1 ∂t v ∂y ( ī) ν+1 ∂x (i) ν ∂y ( ī) ν+1 ∂y (i) ν ∂y ( ī) ν+1 ∂y ( ī) ν                   (45)
for ν = k, k+1, . . . , k+2s+Σ s m=1 Σ 4 j=1 k mj -1, and all the Jacobian matrix components can be computed through Eq. ( 41 The eigenvalues are computed by

|DP -λI| = 0. ( 47
)
Because DP is a 4 × 4 matrix, there are four eigenvalues. If four eigenvalues lie inside the unit circle, then the period-1 motion is stable. If one of them lies outside the unit circle, the periodic motion is unstable. Namely, the stable, periodic motion requires the eigenvalues to be

|λ j | < 1 (j = 1, 2, 3, 4). ( 48
)
If the magnitude of one of the eigenvalues is greater than one,

|λ j | > 1 (j ∈ {1, 2, 3, 4}) (49)
the periodic motion is unstable. For |λ j | < 1 (j = 3, 4) and real λ j (j = 1, 2),

(i) if max{λ j , j = 1, 2} = 1, min{λ j , j = 1, 2} ∈ (-1, 1), ( 50
)
then the saddle-node (SN) bifurcation occurs;

(ii) if min{λ j , j = 1, 2} = -1, max{λ j , j = 1, 2} ∈ (-1, 1), ( 51
)
then the period-doubling bifurcation occurs. For |λ j | < 1 (j = 3, 4) and complex λ j (j = 1, 2),

|λ j | = 1 (j = 1, 2), (52)
then the Neimark bifurcation occurs.

Using the mapping structure in Eq. ( 32), all the periodic motions for the entire range of excitation frequency can be determined analytically by the corresponding governing equations similar to Eqs. ( 42) and ( 43). The mapping structure gives a set of nonlinear algebraic equations, which Table 2. The summary of excitation frequency for impacting chatter (m 1 = 2, m 2 = 1, r 1 = r 2 = 0.6, k 1 = 30, k 2 = 20, Q 0 = 35, d = 1.0 and e = 0.7). can be solved by the Newton-Raphson method.

Mapping

To make the numerical computation convergent fast, the initial guess solutions should be chosen in the appropriate convergent domain. Once the first solution is obtained, the remaining solutions varying with parameters can be determined through the corresponding mapping structure. For the analytical prediction of periodic motions, the parameters ( (2)

m 1 = 2, m 2 = 1, r 1 = r 2 = 0.6, k 1 = 30, k 2 = 20, Q 0 = 30, d = 1.0, e = 0.
-1.5 -1.0 1.0 1.5 SN SN SN AGB SGB k x 5 k x + 1 k x + 4 k x + 2 k x + 3 k x + Excitation Frequency, Ω 2.0 2.2 2.4 2.6 2.8 Switching Velocity, y k (1) -4.0 -2.0 0.0 2.0 4.0 (1) k y (1) 5 k y + (1) 1 k y + (1) 4 k y + (1) 2 k y + (1) 3 k y + SN SN SN AGB SGB (a) (b) Excitation Frequency, Ω 2.0 2.2 2.4 2.6 2.8 Switching Velocity, y k (2) -4.0 -2.0 0.0 2.0 4.0 (2) k y (2) 5 k y + (2) 1 k y + (2) 4 k y + (2) 2 k y + (2) 3 k y + SN SN SN AGB SGB Excitation Frequency, Ω 2.0 2.2 2.4 2.6 2.8 Switch Phase, mod(Ωt k , 2π) 0.0 2.0 4.0 6.0 SN SN SN AGB SGB k t 5 k t + 1 k t + 4 k t + 2 k t + 3 k t + (c) (d)
k = x (1)
k ± d/2. The corresponding switching phase is presented in Fig. 6(d). The switching sets are recorded before impact rather than after impact. The range for the stable, symmetric motion of P 64 2 32 2 is Ω ∈ {(2.457, 2.748), (1, 946, 2.028)}, which is represented by thick curves. The unstable symmetric periodic motion for such a mapping structure lies in Ω ∈ (2.028, 2.457), depicted by the dashed curves.

For the same mapping structure, the asymmetric, periodic motion exists in the range of Ω ∈ (2.380, 2.457), and there are two branches of solutions. The two branches of asymmetric solutions are presented through the dotted and thin curves.

The saddle-node bifurcation of the symmetric, periodic motion occurs at Ω ≈ 2.457 and 2.748 and the saddle-node bifurcation of the asymmetric periodic motion is Ω ≈ 2.457. Such a saddle-node bifurcation causes the symmetry break of impacting chatter motion. The grazing bifurcations for the symmetric and asymmetric periodic motions are Ω ≈ 2.380 and 1.946. For the grazing periodic motion, the eigenvalues show the two periodic motions are stable, i.e. the eigenavalues of the grazing periodic motion lie in the unit circle. However, the grazing causes the mapping structure to be changed. Thus the motion relative to the old mapping structure disappears. In a similar fashion, all the other periodic motion of impacting chatter can be predicted analytically.

Using the same parameters, the analytical prediction of periodic motions with stick is given herein. The corresponding switching sets varying with excitation frequency are presented in Figs. 7(a)-7(d) for the mapping structure of P 654 m 312 m (m = 24, 25, . . . , 28). Because the computational time is very consumed, the impacting chatters with stick are predicted only in the range of excitation frequency Ω ∈ (0.455, 1.39859). The stick bifurcation (SB) occurs at Ω ≈ 1.39859. Once the stick motion appears, the impacting chatter with stick for P 645 26 312 26 are in the range of Ω ∈ (1.328, 1.39859) and (0.635, 0.6789). The stick motion of P 654 27 312 27 lies in two regions of Ω ∈ (1.264, 1.328) and (0.680, 0.9068). The stick motion of P 645 28 312 28 is in Ω ∈ (0.9068, 1.264). The stick motions of P 645 25 312 25 and P 645 24 312 24 are in Ω ∈ (0.612, 0.6789) and (0.455, 0.612), respectively. In stick motion computation, the computational criterion

(i.e. | ẋ(1) k - ẋ(2) k | < 10 -6
) was embedded. If the tolerance is set much smaller, the analytical and numerical predictions cannot be matched sometimes. Such a problem may be caused by the computational error. For a higher accuracy, a better algorithm should be developed. In the other ranges of the stick motion, the analytical prediction can be carried out analytically.

Numerical Simulations

From the analytical prediction, periodic motions for impacting chatter with/without stick will be simulated numerically. All the responses of the first and second gears are represented by the thick and thin curves, respectively. The left and right sides of the second gear are presented by the intermediate thick curves. The switching points before and after impacts are represented by circles. The hollow and gray-filled circles are for the first and second gears, accordingly. The switching points for stick are represented by the large filled circles.

For system parameters (m 1 = 2, m 2 = 1, r 1 = r 2 = 0.6, k 1 = 30, k 2 = 20, Q 0 = 30, d = 1.0 and e = 0.7), consider the excitation frequency Ω = 2.72. The analytical prediction gives the initial condition (i.e. t 0 ≈ 2.0559169, x

(1) 0 68706528 andy (2) 0 ≈ 4.5219068) for the stable periodic impacting chatter of P 64 2 32 2 . The displacement and velocity timehistories, trajectories of the first and second gears in phase plane are presented in Figs. 8(a)-8(d). The displacement responses of two gears are presented in Fig. 8(a). The motion starts at the right side (Rside) of the second gear. The first gear with the second gear has two impacts. Afterwards, the two gears have a free-flying motion. Once the first gear arrives to the left side (L-side) of the second gear, the next two impacts between the two gears occurs. After these two impacts, the two gears have another freeflying motion, and the two gears meet on the right side of the second gear. Such a motion process will repeat and form a periodic motion. Due to impacts between two gears on the left and right sides, the impacting velocities become discontinuous. Thus, the velocity time-histories for two gears are presented in Fig. 8(b). Such discontinuity is observed at the impact locations. After the right side impact, the first gear velocity is less than the second gear velocity (i.e. y the second gear can be obtained by adding the two offsets.

≈ 1.0505790, x (2) 0 ≈ 0.5505790, y (1) 0 ≈ 1.
To understand such impacting chatter in the gear transmission system, the relative force responses should be discussed. After the first rightside impact, to make the repeated impact occur, the relative force should be greater than zero (i.e. g (1)

ky

(2) k < 0). Thus, just after the impact, the relative displacement is less than half a gap (i.e. z

(1) k = x (1) k -x (2) k < d/2).
To allow the first gear to catch up with the second one, the relative force should be greater than zero (i.e. g

(1) α (t k ) > 0) before the next impact. Otherwise, the next impact on the righthand side will not happen. In a similar fashion, just after the first impact on the left-side the relative force should be less than zero (i.e. g

(1) α (t k ) < 0) to achieve the second impact. From Figs. 9(a) and 9(b), such a relation is observed. Based on the perspective of the second gear, the relative force criterion of the periodic motion is described through z (2) = x (2)x (1) = -z (1) , and the relative velocity and force for the second gear should be opposite to the first gear. So the relative force for the second gear is presented in Figs. 9(a) and 9(b) as well.

The impacting chatter with stick in the gear transmission system is very interesting. Such a periodic motion will be discussed through the mapping structure of P 5 25 312 25 64 . The impacting chatter with stick relative to P 5 25 312 25 64 is equivalent to the one relative to P 645 25 312 25 . Just the starting points are from different switching points. As before, the same parameters will be adopted for this illustration with Ω = 0.5. From the analytical prediction, the initial condition is selected from the switching point t 0 ≈ 4.8888960, x

(1) 0 ≈ -0.7387947, x

(2) 0 ≈ -0.2387947, and y

(1) 0 = y (2) 0 ≈ -0.2754161. From such an initial condition, the stick motion will appear on the left side of the second gear. After the stick vanishes, the first gear moves to the right side of the second gear through the free-flying motion. Because of different velocities of the two gears, the impacting chatter occurs before the stick appears. For this motion, once y

(1)

k = y (2)
k , the stick motion will be observed. In Fig. 10, the displacement, velocity and trajectories in phase planes are presented. The shaded area is used for stick motion, and the Letters "L" and "R" represent the left-side and right-side of the second gear, respectively. In Fig. 10(a), the displacement responses of the two gears are presented. The stick motion on the left and right sides of the second gear are observed. The impacting chatter with 26 impacts on each side is observed. To further confirm the impacts, the velocity responses for the two gears are presented in Fig. 10(b). For further illustration of the impacting chatter with stick, the trajectories of the first and second gears in phase plane are plotted in Figs. 10(c) and 10(d), respectively.

For a better understanding of impacting chatter with stick, the relative force and jerk for the two gears are illustrated in Fig. 11 β (t k ), α, β ∈ {1, 2, 3} and α = β). From [START_REF] Luo | Nonlinear dynamics of a gear transmission dynamical system, Part I: Mechanism of impacting chatter and stick[END_REF], the onset of the stick motion requires the two forces g α (t k ) = 0, which is observed in Figs. 11(a) and 11(b). To correspond to the trajectory in phase plane, the relative force distribution along the displacement of the first gear is shown in Fig. 11(c). The aforementioned conditions for the onset and vanishing of the stick motion are observed. However, this condition of g (1) α (t) < 0 for t > t k + ε, the corresponding relative velocity and acceleration are greater than zero, which leads to the relative displacement satisfying z (1) = x (1)x (2) < d/2. In other words, the two gears lie in the state of free-flying motion, and the stick disappears. The stick motion disappearance on the left side can be discussed in a similar manner. The jerk condition of J

(1) α (t k ) > 0 leads to z (1) = x (1)x (2) > -d/2. Therefore, the two gears start to be separated from the left side of the second gear, and the stick motion disappears.

Conclusions

In this paper, the analytical prediction of periodic motions relative to the impacting chatter with and without stick in the gear transmission system is completed through the mapping structure. The corresponding local stability and bifurcation analysis are carried out, and the grazing and stick conditions are used for the existence of periodic motions. Numerical simulations are performed to illustrate the periodic motions and stick motion criteria. The prediction of chaotic motions in such a gear system will be discussed in another paper.

Fig. 1 .

 1 Fig. 1. Basic mappings for impacting chatter in (a) the absolute frame and (b) the relative frame. The straight line with arrow represents an impact on the boundary.

  . The mappings relative to the stick switching planes only are sketched in Fig. 2. In Fig. 2(a), only two stick mappings (P 1 and P 4 ) are new, and the other four mappings are the same as in Fig. 1. In the relative frame, the switching planes are points. The stick mapping is very difficult to be illustrated. The other mappings are presented in Figs. 2(b). The mappings based on the sticking and impacting switching planes are presented in Figs. 3(a) and 3(b).

Fig. 2 .

 2 Fig. 2. Mappings from stick switching planes: (a) absolute motion and (b) relative motion.

Fig. 3 .

 3 Fig. 3. Mappings between switching planes for stick and impacting.

Fig. 4 .

 4 Fig. 4. Mapping structures for (a) impacting chatter and (b) impacting chatter with stick motion of two gear systems.

Fig. 5 .

 5 Fig. 5. Bifurcation scenario: (a) switching displacement of the first gear, (b) switching velocity of the first gear, (c) switching velocity of the second gear, and (d) switching phase (m 1 = 2, m 2 = 1, r 1 = r 2 = 0.6, k 1 = 30, k 2 = 20, Q 0 = 30, d = 1.0 and e = 0.7).

Fig. 6 .

 6 Fig. 6. The analytical prediction of a periodic motion of impacting chatter P 64 2 32 2 : (a) switching displacement and (b) switching velocity of the first gear, (c) switching velocity of the second gear, and (d) switching phase (m 1 = 2, m 2 = 1, r 1 = r 2 = 0.6, k 1 = 30, k 2 = 20, Q 0 = 30, d = 1.0 and e = 0.7).

Fig. 7 .

 7 Fig. 7. The analytical prediction of a periodic motion of impacting chatter with stick P 645 m 312 m (m = 24, 25, . . . , 28): (a) switching displacement and (b) switching velocity of the first gear, (c) switching velocity of the second gear, and (d) switching phase (m 1 = 2, m 2 = 1, r 1 = r 2 = 0.6, k 1 = 30, k 2 = 20, Q 0 = 30, d = 1.0 and e = 0.7).

kFig. 8 .

 8 Fig. 8. The periodic motion of impacting chatter P 65 2 32 2 : (a) displacement time-history and (b) velocity time-history, (c) trajectories of the first gear, (d) trajectories of the second gear in phase plane (m 1 = 2, m 2 = 1, r 1 = r 2 = 0.6, k 1 = 30, k 2 = 20, Q 0 = 30, d = 1.0, e = 0.7 and Ω = 2.72). The initial condition is t 0 ≈ 2.0559169, x (1) 0 ≈ 1.0505790, x (2) 0 ≈ 0.5505790, y (1) 0 ≈ 1.68706528 and y (2) 0 ≈ 4.5219068.

α

  (t k ) > 0). Since the absolute velocity relation (y (1) k < y (2) k ) holds just after impact, it means the relative velocity is less than zero ( ż(1) k = y

Fig. 9 .

 9 Fig. 9. Relative forces for the periodic motion of impacting chatter P 65 2 32 2 : (a) time-history and (b) distribution along the displacement of the first gear (m 1 = 2, m 2 = 1, r 1 = r 2 = 0.6, k 1 = 30, k 2 = 20, Q 0 = 30, d = 1.0, e = 0.7 and Ω = 2.72). The initial condition is t 0 ≈ 2.0559169, x (1) 0 ≈ 1.0505790, x (2) 0 ≈ 0.5505790, y (1) 0 ≈ 1.68706528 and y (2) 0 ≈ 4.5219068.

Fig. 10 .

 10 Fig. 10. The periodic motion of impacting chatter with stick P 645 25 312 25 : (a) displacement time-history and (b) velocity time-history, (c) trajectories of the first gear, (d) trajectories of the second gear in phase plane (m 1 = 2, m 2 = 1, r 1 = r 2 = 0.6, k 1 = 30, k 2 = 20, Q 0 = 30, d = 1.0, e = 0.7 and Ω = 0.5). The initial condition is t 0 ≈ 4.8888960, x (1) 0 ≈ -0.7387947, x (2) 0 ≈ -0.2387947 and y (1) 0 = y (2) 0 ≈ -0.2754161.

  k ) < 0) on the right-side for α = 1, 2 and g

α

  (t k ) < 0 (or g

αFig. 11 .

 11 Fig. 11. Forces for the periodic motion of impacting chatter with stick P 645 25 312 25 : (a) relative force time-history and (b) internal forces for stick, (c) distributions along the displacement and (d) the relative jerk of the first gear (m 1 = 2, m 2 = 1, r 1 = r 2 = 0.6, k 1 = 30, k 2 = 20, Q 0 = 30, d = 1.0, e = 0.7 and Ω = 0.5). The initial condition is t 0 ≈ 4.8888960, x (1) 0 ≈ -0.7387947, x (2) 0 ≈ -0.2387947 and y (1) 0 = y (2) 0 ≈ -0.2754161.

  g

  (i) α (t k ) and g

  (i) β (t k ) to have the same sign. To obtain the stick motion in Ω 1 requires g

  k ) < 0 and g

  k ) < 0; and to obtain the stick motion in Ω 3 requires g

  k ) > 0. Such a relative force characteristic is presented in Fig.11(a). InFig. 11(b), the internal forces in the stick domains for two gears are presented because the internal forces between the two gears only exist for the stick. When the stick motion disappears, the internal forces vanish. In other words, the internal forces are zero at the vanishing of the stick motion. The vanishing condition for stick is g

  ) = 0 is a necessary condition. The sufficient condition for the vanishing of the stick motion requires J(1) α (t k ) < 0 (or J

  k ) > 0 (or J

  k ) < 0) on the left side. Therefore, in Fig.11(d), the jerk J

  the first gear is presented, and such a condition is observed. For instance, when the stick vanishes at z k = d/2, the relative velocity and acceleration are zero (i.e. z

  ). The variational equation for a set of switching points {y *

			k+2s+Σ s m=1 Σ 4 j=1 k mj -1 }
	is
	∆y k+2s+	P s m=1	P 4 j=1 k mj = DP(y * k )∆y k	(46)

k , . . . , y *

Table 1 .

 1 The summary of excitation frequency for impacting chatter with stick (m 1 = 2, m 2 = 1, r 1 = r 2 = 0.6, k 1 = 30, k 2 = 20, Q 0 = 35, d = 1.0 and e = 0.7).

	Mapping Structure	Excitation Frequency
	P 62 26 135 26 4	P (T )	( 1 .382, 1.397)
	P 62 27 135 27 4	P (T )	( 1 .264, 1.382)
	P 62 28 135 28 4	P (T )	( 0 .908, 1.264)
	P 62 27 135 27 4	P (T )	( 0 .680, 0.908)
	P 62 26 135 26 4	P (T )	( 0 .635, 0.680)
	P 62 25 135 25 4	P (T )	( 0 .613, 0.635)
	P 62 24 135 24 4	P (T )	( 0 .455, 0.613)
	P 62 25 135 25 4	P (T )	( 0 .346, 0.455)
	P 62 24 135 24 4	P (T )	( 0 .317, 0.346)
	P 62 23 135 23 4	P (T )	( 0 .216, 0.317)
	P 62 k 135 k 4	P (T )	( 0 .000, 0.216)

  7) are used. The analytical prediction of the periodic motion of P 64 2 32 2 gives switching sets varying with excitation frequency, as shown in Fig.6. Because of the page limitation, the eigenvalues varying with

	(1)
	k
	Switching Displacement x
	2.0	2.2	2.4	2.6	2.8
		Excitation Frequency, Ω

Appendix

With an initial condition (t, x (i) , ẋ(i) ) = (t k , x

k ), the general solutions for two gear oscillators for three regions are given as follows:

where

where

where

where

where .11)


Alp Ár

Rich Árd

M Ész

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

ON THE VARIATIONAL FORMULATION OF SOME STATIONARY SECOND ORDER MEAN FIELD GAMES SYSTEMS

Keywords: AMS Subject Classification: 49K20, 35J47, 35J50, 49N70, 91A13 Mean Field Games, weak solutions, variational formulation, density constraints, general couplings

. Furthermore, our approach can be used to obtain solutions of systems with multiple populations.

Introduction

In this article we consider the stationary Mean Field Game (MFG) system

(MFG 1 ) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ -∆u + H(⋅, ∇u) + λ = f (x, m) in Ω, ∇u ⋅ n = 0 on ∂Ω,
-∆mdiv m∇ ξ H(⋅, ∇u) = 0, in Ω,

(∇m + m∇ ξ H(⋅, ∇u)) ⋅ n = 0, on ∂Ω, Ω m dx = 1, m(x) > 0 in Ω,
introduced by Lasry and Lions in [START_REF] Lasry | Jeux à champ moyen I. Le cas stationnaire[END_REF][START_REF] Lasry | Mean field games[END_REF], which models the equilibrium configuration of an ergodic stochastic symmetric differential game with a continuum of small, indistinguishable players (see e.g. [START_REF] Feleqi | The derivation of ergodic mean field game equations for several populations of players[END_REF]). In the above system, Ω ⊆ R d is a bounded domain with a smooth boundary, n is the outward normal to ∂Ω, f ∶ Ω × [0, +∞[→ R is the so-called local coupling function, Ω × R d ∋ (x, ξ) ↦ H(x, ξ) ∈ R is the Hamiltonian, and the variables m ∶ Ω → [0, +∞), u ∶ Ω → R and λ ∈ R represent the stationary equilibrium configuration of the players, the equilibrium cost of a typical player, and the ergodic constant, respectively. We have taken Neumann boundary conditions but our results admit natural versions for the more standard case of periodic boundary conditions. An interesting feature of the solutions of (MFG 1 ) is their connection with the long time behaviour of the solutions of time-dependent MFGs. We refer the reader to [START_REF] Cardaliaguet | Long time average of mean field games[END_REF][START_REF] Cardaliaguet | Long time average of mean field games with a nonlocal coupling[END_REF] for some results justifying rigorously this relation in some special cases. The numerical resolution of stationary second order MFGs has been studied in [START_REF] Achdou | Mean field games: Numerical methods[END_REF][START_REF] Cacace | A generalized newton method for homogenization of hamilton-jacobi equations[END_REF][START_REF] Almulla | Two numerical approaches to stationary mean-field games[END_REF][START_REF] Briceno-Arias | Proximal methods for mean field games with local couplings[END_REF].

Existence and uniqueness results for system (MFG 1 ) have been investigated by several researchers using Partial Differential Equation (PDE) techniques, starting with the first papers [START_REF] Lasry | Jeux à champ moyen I. Le cas stationnaire[END_REF][START_REF] Lasry | Mean field games[END_REF] in the framework of weak solutions. The reader is referred to [START_REF] Cirant | Multi-population mean field games systems with Neumann boundary conditions[END_REF][START_REF] Ferreira | Existence of weak solutions to stationary mean-field games through variational inequalities[END_REF][START_REF] Bardi | Nonlinear elliptic systems and mean-field games[END_REF][START_REF] Cirant | Stationary focusing mean-field games[END_REF] for other subsequent results on the existence of weak solutions. In addition, under different assumptions on the coupling function f and the growth of H, the existence, and uniqueness, of smooth solutions have been analysed in [START_REF] Gomes | On the existence of classical solutions for stationary extended mean field games[END_REF][START_REF] Gomes | Existence for stationary mean-field games with congestion and quadratic Hamiltonians[END_REF][START_REF] Pimentel | Regularity for second order stationary mean-field games[END_REF]. See also [START_REF] Gomes | A-priori estimates for stationary mean-field games[END_REF], where several a priori estimates for smooth solutions of stationary second order MFGs are established.

In this article we focus our attention on the proof of the existence of weak solutions of (MFG 1 ) by variational techniques. Indeed, as pointed out already in [START_REF] Lasry | Mean field games[END_REF], system (MFG 1 ) can be seen, formally, as the first order optimality condition of an associated variational problem, involving a PDE constraint for the variable m. It turns out that u and λ in (MFG 1 ) correspond to the Lagrange multipliers associated to the PDE constraint for m and the condition ∫ Ω m(x) dx = 1, respectively. Given q > d ≥ 2, where d is the space dimension, and setting q ′ ∶= q (q -1), we prove the above assertion for Hamiltonians H growing as ⋅ q ′ . Even if this growth condition is restrictive, but crucial for our arguments, the main interest of this variational technique is that it allows to prove the existence of weak solutions of (MFG 1 ) for a rather general class of coupling functions f in a straightforward manner. Indeed, as we will show in Section 3, f does not need to be monotone (see also [START_REF] Cirant | Stationary focusing mean-field games[END_REF][START_REF] Cirant | Singular mean-filed games[END_REF] for some recent results in this direction) and, moreover, we can prove the existence of solutions of variations of system (MFG 1 ) involving couplings which can also depend on the distributional derivatives of m. As a matter of fact, our results are valid, for terms in the r.h.s. of the first equation in (MFG 1 ) which can be identified with the derivative of a function F ∶ W 1,q (Ω) → R which is Gâteaux dfferentiable and weakly lower-semicontinuous.

Our approach follows closely the one in [START_REF] Mészáros | A variational approach to second order mean field games with density constraints: the stationary case[END_REF], which considers in addition a density constraint in order to model strong congestion effects (see [START_REF] Santambrogio | A modest proposal for MFG with density constraints[END_REF][START_REF] Cardaliaguet | First order mean field games with density constraints: pressure equals price[END_REF]). In that article, the existence of solutions (m, u, p, µ, λ) ∈ W 1,q (Ω) × W 1,q ′ (Ω) × M (Ω) × M (Ω) × R, where M (Ω) is the set of Radon measures on Ω, to the system

(1.1) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ -∆u + 1 q ′ ∇u q ′ + µ -p + λ = f (x, m) in Ω, ∇u ⋅ n = 0 on ∂Ω, -∆m -div ∇u q ′ -2 ∇um = 0, in Ω, ∇m ⋅ n = 0, on ∂Ω, Ω m dx = 1, 0 ≤ m(x) ≤ 1 in Ω, spt(p) ⊆ {x ∈ Ω ; m(x) = 1}, p ≥ 0 in Ω, spt(µ) ⊆ {x ∈ Ω ; m(x) = 0}, µ ≥ 0 in Ω, is established if f (x, ⋅) is non-decreasing. If q > d
, the result is proved with a variational approach. On the other hand, when 1 < q ≤ d a penalization argument allows to prove the existence also in this case. In the present article, when q > d, we improve the results in [START_REF] Mészáros | A variational approach to second order mean field games with density constraints: the stationary case[END_REF], and we show the existence of solutions of

(MFG 2 ) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ -∆u + H(⋅, ∇u) -p + λ = f (x, m) in Ω, ∇u ⋅ n = 0 on ∂Ω, -∆m -div m∇ ξ H(⋅, ∇u) = 0, in Ω, (∇m + m∇ ξ H(⋅, ∇u)) ⋅ n = 0, on ∂Ω, Ω m dx = 1, 0 < m(x) ≤ κ(x) in Ω, spt(p) ⊆ {x ∈ Ω ; m(x) = κ(x)}, p ≥ 0 in Ω.
Note that in (MFG 2 ) a more general Hamiltonian is considered and the density constraint m ≤ 1 is replaced by m ≤ κ, where κ ∈ W 1,q (Ω). Most importantly, f does not need to be monotone and, using the Harnack's inequality proved in [START_REF] Trudinger | Linear elliptic operators with measurable coefficients[END_REF] (see also [START_REF] Bogachev | Fokker-Planck-Kolmogorov equations[END_REF]) for elliptic equations in divergence form, we show that the density m is strictly positive, which implies that µ in (1.1) is identically zero. Using the existence of solutions of the variational problem associated to (MFG 2 ), which can be proved easily, we prove that the variational problem associated to (MFG 1 ) admits at least one solution. This crucial fact is the key to show the existence of solutions of (MFG 1 ). The paper is organized as follows. In Section 2 we begin with some preliminaries which allow us to characterize the subdifferential of the cost functionals appearing in the optimization problems associated to (MFG 1 )-(MFG 2 ). This analysis extends the one in [MS15, Section 2]. Section 3 is the core of the article. We prove the existence of solutions of the variational problems associated to (MFG 1 )-(MFG 2 ) and we establish the corresponding optimality conditions, which provide the existence of solutions of (MFG 1 )-(MFG 2 ). We present a detailed discussion concerning the generality of the coupling term, which, as we have explained before, is the main feature of this approach. We also prove, by a bootstrapping argument, additional regularity for the weak solutions. In Section 4, we present some simple applications of our results to the study of multi-populations MFG systems (see e.g. [START_REF] Cirant | Multi-population mean field games systems with Neumann boundary conditions[END_REF][START_REF] Bardi | Nonlinear elliptic systems and mean-field games[END_REF]). Finally, in the appendix, we prove the strict positivity of the densities m appearing in (MFG 1 )-(MFG 2 ) as a consequence of the Harnack's inequality in [START_REF] Trudinger | Linear elliptic operators with measurable coefficients[END_REF] and the assumed regularity of the boundary ∂Ω.

Preliminary results

In the entire article, we will assume that Ω ⊆ R d (d ≥ 2) is a non-empty, bounded open set with a C 1,1 boundary ∂Ω. This regularity assumption is equivalent to a uniform interior and exterior ball condition (see for instance [Dal14, Theorems 1.8-1.9]) and allows us to use the classical Sobolev inequalities. The vector n will denote the outward normal to ∂Ω. Given r ∈ [1, +∞] and ℓ ∈ N we will denote by ⋅ r and ⋅ ℓ,r the standard norms in L r (Ω) (or in L r (Ω) d ) and W ℓ,r (Ω), respectively.

Let q > d ≥ 2. Our aim in this section is to provide a characterization of the subdifferential of the convex functional

B q ∶ W 1,q (Ω) × L q (Ω) d → R ∶= R ∪ {+∞}, defined as (2.1) B q (m, w) ∶= Ω b q (x, m(x), w(x)) dx, where b q ∶ Ω × R × R d → R is given by (2.2) b q (x, m, w) ∶= ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ mH * (x, -w m) if m > 0, 0 if (m, w) = (0, 0), +∞ otherwise.
In view of the assumptions below, the function B q is well defined and convex (see Remark 2.1 and Theorem 2.4). It will appear in the cost functional of an optimization problem whose first order optimality condition has the form of an MFG system. In (2.2), for every x ∈ Ω the function

H * (x, ⋅) ∶ R d → R is the Legendre-Fenchel transform of H(x, ⋅), where H ∶ Ω × R d → R is a
continuous function, that we will call the Hamiltonian, which is assumed to be strictly convex and differentiable in its second variable and satisfies a polynomial growth condition in terms of q ′ ∶= q (q -1): there exist

C 1 , C 2 > 0 such that (2.3) 1 q ′ C 1 ξ q ′ -C 2 ≤ H(x, ξ) ≤ C 1 q ′ ξ q ′ + C 2 , ∀ x ∈ Ω, ξ ∈ R d .
Using the definition of H * (x, ⋅), an easy computation shows that (2.3) implies

(2.4)

C 1-q 1 q η q -C 2 ≤ H * (x, η) ≤ C q-1 1 q η q + C 2 , ∀ x ∈ Ω, η ∈ R d .
Note that since H(x, ⋅) is strictly convex and differentiable, we have that H * (x, ⋅) is also strictly convex and differentiable (see e.g. [START_REF] Rockafellar | Convex analysis[END_REF]Theorem 26.3]). We denote their gradients by ∇ ξ H(x, ⋅) and ∇ η H * (x, ⋅), respectively. Moreover, classical results in convex analysis show that for any x ∈ Ω (2.5)

∇ ξ H(x, ⋅) -1 (η) = ∇ η H * (x, η) ∀ η ∈ R d .
We now prove an elementary result which will be useful later. In the remainder of this article will denote by C > 0 a generic constant which can change from line to line.

Lemma 2.1. Under (2.3) we have that H * (⋅, ⋅), ∇ ξ H(⋅, ⋅) and ∇ η H * (⋅, ⋅) are continuous. Moreover, if

β ∈ L q ′ (Ω) d we have that ∇ ξ H(⋅, β(⋅)) ∈ L q (Ω) d . Analogously, if v ∈ L q (Ω) d then ∇ η H * (⋅, v(⋅)) ∈ L q ′ (Ω) d . Proof. Let (x n , η n ) and (x, η) in Ω × R d be such that (x n , η n ) → (x, η) as n → ∞. Set ξ n ∶= ∇ η H * (x n , η n ). By definition of H * H * (x n , η n ) = η n ⋅ ξ n -H(x n , ξ n ) ≥ η n ⋅ ξ -H(x n , ξ) ∀ ξ ∈ R d .
Since η n is convergent, hence bounded, the first inequality in (2.3) shows that ξ n is bounded. Let ξ be a limit point of ξ n . Then, using the continuity of H and passing to the limit, up to some subsequence, we get that

η ⋅ ξ -H(x, ξ) ≥ η ⋅ ξ -H(x, ξ) ∀ ξ ∈ R d ,
which shows that ξ = ∇H * (x, η) and any limit point of H * (x n , η n ) is equal to H * (x, η). The continuity of ∇H(⋅, ⋅) follows by the symmetric argument. Finally, if

β ∈ L q ′ (Ω) d , setting x ∈ Ω ↦ η(x) ∶= ∇H(x, β(x))
, which is a measurable function since ∇H(⋅, ⋅) is continuous, we get by convexity that H * (x, η(x)) ≤ β(x) ⋅ η(x) + H * (x, 0), and so, by (2.4), we obtain the existence of C > 0 such that η(x) q ≤ C(β(x) ⋅ η(x) + 1). Using Young's inequality, we get the existence of C > 0 such that η(x) q ≤ C( β(x) q ′ + 1) and thus, integrating in Ω, we obtain that η ∈ L q (Ω). The last assertion follows from an analogous argument.

Regarding the dependence of H on the space variable x, we will assume that there exists a modulus of continuity which is uniform w.r.t. the second variable, i.e. ∃ ω ∶ [0, +∞) → [0, +∞) such that ω(0) = 0, ω is continuous, non-decreasing and

(2.6) H(x, ξ) -H(y, ξ) ≤ ω( x -y )( ξ q ′ + 1), ∀ x, y ∈ Ω, ξ ∈ R d .
Using that Ω is a compact set, a natural example of a Hamiltonian H satisfying (2.3) and (2.6) is given by H(x, ξ) ∶= b(x) ξ q ′ + c(x) where b, c ∈ C(Ω) and b > 0. Following the analysis in [MS15, Lemma 2.1 and Theorem 2.2], presented in a more particular setting, we shall characterize the subdifferential of B q , defined in (2.1). Recall that given a normed space (X, ⋅ ) and a l.s.c. convex proper function g ∶ X → R, the subdifferential ∂g(x) of f at the point x, consist in the set of all x * ∈ X * such that

g(x) + ⟨x * , y -x⟩ X * ,X ≤ g(y) ∀ y ∈ X.
For the sake of completeness, in order to identify ∂B q , we first state some simple properties of the function b q (x, ⋅, ⋅). Given x ∈ Ω consider the set (2.7)

A q ′ (x) ∶= {(α, β) ∈ R × R d ∶ α + H(x, -β) ≤ 0}.
Since H is continuous and convex w.r.t. its second variable, we have that A q ′ (x) is closed and convex for any x ∈ Ω. Given a subset D of an euclidean space, we denote by χ D its characteristic function (in the sense of convex analysis), i.e. χ D (y) = 0 if y ∈ D and χ D (y) = +∞ otherwise.

Lemma 2.2. For all x ∈ Ω, the function b q (x, ⋅, ⋅) is convex, proper and l.s.c. Its Legendre-Fenchel conjugate and its subdifferential are given by

(2.8) b * q (x, ⋅, ⋅) = χ A q ′ (x) (⋅, ⋅), ∂ (m,w) b q (x, m, w) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ (-H(x, -β x ), β x ) if m > 0, (α, β) ∈ A q ′ (x) if (m, w) = (0, 0), ∅ otherwise,
where, if m > 0, β x ∶= -∇H * (x, -w m).

Proof. Using (2.2), it is straightforward to check that for all x ∈ Ω we have

(2.9) b q (x, m, w) = sup β∈R d {β ⋅ w -mH(x, -β)} ,
and so b q (x, ⋅, ⋅) is convex, proper and l.s.c. Using that b q (x, ⋅, ⋅) is proper, we have that

∂ (m,w) b q (x, m, w) = ∅ if m < 0 or m = 0 and w ≠ 0.
On the other hand, from (2.9), if m ≥ 0, for any x ∈ Ω we have the identity

(2.10) b q (x, m, w) = sup {αm + β ⋅ w ∶ (α, β) ∈ A q ′ (x)} ,
which is checked to hold also when m < 0. Thus, since

A q ′ (x) is closed and convex, b * q (x, ⋅, ⋅) = χ A q ′ (x) (⋅, ⋅). This expression directly yields that ∂ (m,w) b q (x, 0, 0) = A q ′ (x). If m > 0, then b q (x, ⋅, ⋅
) is differentiable and, by a simple computation, we get the expression of its gradient with respect to (m, w).

Remark 2.1. Notice that the equality in (2.10) shows that (x, m, w) ↦ b q (x, m, w) is lowersemicontinuous and so, by [RW98, Example 14.31], we have that b q is a normal integrand. This shows that Ω ∋ x ↦ b q (x, m(x), w(x)) is a measurable function if m and w are measurable (see [RW98, Proposition 14.28]). In particular, the functional B q is well defined.

Let us define

(2.11) Lemma 2.3. The closure of A q ′ in (W 1,q (Ω)) * × L q ′ (Ω) d is given by (2.12)

A q ′ ∶= (a, b) ∈ L ∞ (Ω) × L ∞ (Ω) d ∶ (a(x), b(x)) ∈ A q ′ (x)
A q ′ ∶= (α, β) ∈ M (Ω) × L q ′ (Ω) d ∶ α + H(⋅, -β) ∈ M -(Ω) ,
or equivalently,

(2.13) A q ′ ∶= (α, β) ∈ M (Ω) × L q ′ (Ω) d ∶ α ac + H(⋅, -β) ≤ 0, a.e. in Ω and α s ∈ M -(Ω) ,
where dα = α ac dx + dα s is the Lebesgue decomposition of the measure α w.r.t. the Lebesgue measure restricted to Ω.

Proof. Let us take (α n , β n ) ∈ A q ′ converging to (α, β) in (W 1,q (Ω)) * × L q ′ (Ω). By definition of A q ′ one has Ω α n (x)φ(x) dx ≤ - Ω H(x, -β n (x))φ(x) dx, ∀φ ∈ W 1,q (Ω), φ ≥ 0. Since β n → β in L q ′ (Ω) d , up to some subsequence β n (x) → β(x) a.e. in Ω. By (2.3) we can use Fatou's lemma to obtain that ⟨α, φ⟩ (W 1,q ) * ,W 1,q ≤ - Ω H(x, -β(x))φ(x) dx, ∀φ ∈ W 1,q (Ω), φ ≥ 0.
Using (2.3) again, we obtain that α -C 2 defines a nonpositive distribution, hence by [Sch66, Théorème V] α can be identified with an element of M (Ω).

Conversely let us take (α, β) belonging to the right-hand-side (r.h.s.) of (2.12), or equivalently to the r.h.s. of (2.13). Analogously to [START_REF] Mészáros | A variational approach to second order mean field games with density constraints: the stationary case[END_REF], we construct different approximations for α ac and β on the one hand and for α s on the other hand. For R > 0 and

x ∈ R d we set B R (x) = {y ∈ R d ∶ y -x < R}. Consider a mollifier η ∶ R d → R satisfying that η ∈ C ∞ c (R d ), η ≥ 0, R d η(x) dx = 1, spt(η) ⊆ B 1 (0) and η(x) = η(-x) for all x ∈ R d . Now, for ε > 0 set Ω ε ∶= {x ∈ Ω ∶ dist(x, ∂Ω) > ε}, η ε (x) ∶= 1 ε d η(x ε)
, and for all x ∈ Ω and i = 1, ..., d, let us define

αε (x) ∶= Ω η ε (x -y)α ac (y) dy½ Ωε (x), βi ε (x) ∶= Ω η ε (x -y)β i (y) dy½ Ωε (x).
As ε ↓ 0, we have that αε

→ α ac in L 1 (Ω), βε → β in L q ′ (Ω) d and (α ε , βε ) ∈ L ∞ (Ω) × L ∞ (Ω) d .
Multiplying the inequality in (2.13) by η ε , integrating and using Jensen's inequality yield

(2.14) αε (x) + H(x, -βε (x)) ≤ δε (x), ∀x ∈ Ω, where δε (x) ∶= Bε(x) [H(x, -β(y)) -H(y, -β(y))] η ε (x -y) dy½ Ωε (x).
Note that (2.6) implies that (2.15) δε (x) ≤ ω(ε)

Bε(x) 1 + β(y) q ′ η ε (x -y) dy½ Ωε (x) ≤ ω(ε) ε d η ∞ 1 + β q ′ 1 ,
and so δε

∈ L ∞ (Ω). Using that g ε (⋅) ∶= ∫ Ω 1 + β(y) q ′ η ε (⋅ -y) dy½ Ωε (⋅) converges in L 1 (Ω) to 1 + β(⋅) q ′ ∈ L 1 (Ω)
, extracting a subsequence, the first inequality in (2.15) implies that δε (⋅) → 0 a.e. in Ω. Since for ε ∈ (0, 1) we have that δε

≤ ω(1)g ε , we get from [EG92, Chapter 1.3, Theorem 4] that δε → 0 in L 1 (Ω). This implies that (α ε , βε ) ∶= (α ε -δε , βε ) ∈ A q ′ and (α ε , βε ) → (α ac , β) in L 1 (Ω) × L q ′ (Ω) d .
Now, in order to approximate the singular part α s , for x ∈ Ω and ε > 0 let us define ρ x ε ∶=

½ Bε (x)∩Ω B ε (x) ∩ Ω and αs ε (y) ∶= Ω ρ x ε (y) dα s (x) ∀ y ∈ Ω,
which is a non-positive function. Arguing exactly as in the proof of [MS15, Lemma 2.1] we get that the uniform interior ball assumption on the boundary ∂Ω implies that αs

ε ∈ L ∞ (Ω). Using that ρ x ε → δ x in M (Ω), as ε ↓ 0, it is straightforward to show that αs ε → α s in M (Ω). Therefore, the sequence (α s ε + αε , βε ) → (α, β) in M (Ω) × L q ′ (Ω) d and (α s ε + αε , βε ) ∈ A q ′ for
all ε > 0. Using that q > d, and so W 1,q (Ω) ↪ C(Ω) by the Sobolev embedding, we have that

(α s ε + αε , βε ) → (α, β) weakly in (W 1,q (Ω)) * × L q ′ (Ω) d .
Since A q is convex, its closure w.r.t. the weak and strong topologies coincide. The result follows.

For a given representative of m

∈ W 1,q (Ω) in C(Ω), we denote {m = 0} ∶= x ∈ Ω ∶ m(x) = 0 and {m > 0} ∶= x ∈ Ω ∶ m(x) > 0 .
Theorem 2.4. The following assertions hold true: (i) The functional B q , defined in (2.1), is convex, l.s.c. and

B * r (α, β) = χ A r ′ (α, β) for all (α, β) ∈ (W 1,q (Ω)) * × L q ′ (Ω) d . (ii) Let (m, w) ∈ W 1,q (Ω) × L q (Ω) d and suppose that B q (m, w) < ∞. Then, if v ∶= (w m)½ {m>0} ∉ L q (Ω) d we have that ∂B q (m, w) = ∅.
Otherwise, B q is subdifferentiable at (m, w) and

(2.16) ∂B q (m, w) = (α, β) ∈ A q ′ ∶ α {m > 0} = -H(⋅, ∇H * (⋅, -v)) and β {m > 0} = -∇H * (⋅, -v) .
In particular, the singular part of α in (2.16), w.r.t. to the Lebesgue measure, is concentrated in {m = 0}.

Proof. Since the arguments are similar to those in the proof of [MS15, Theorem 2.2], we only sketch the main ideas. First, truncating the sets A q ′ (x), defined in (2.7), by setting for k ∈ N

A q ′ ,k (x) ∶= (a, b) ∈ A q ′ (x) ; a ≥ -k, max i=1,⋯,d b i ≤ k ,
using Lemma 2.2 and the monotone convergence theorem, we have that

B q (m, w) = Ω sup (a,b)∈A q ′ (x) {am(x) + b ⋅ w(x)} dx = lim k→∞ Ω sup (a,b)∈A q ′ ,k (x) {am(x) + b ⋅ w(x)} dx.
Characterizing the point-wise optimizers (α(x), β(x)) ∈ A q ′ ,k (x) in the last expression, it easy to see that Ω ∋ x ↦ (α(x), β(x)) is a measurable function, which by definition belongs to

L ∞ (Ω) × L ∞ (Ω) d .
Using this fact and the monotone convergence theorem once again, we find that B q (m, w) = sup

(α,β)∈A q ′ Ω {α(x)m(x) + β(x) ⋅ w(x)} dx = sup (α,β)∈A q ′ Ω m(x) dα(x) + Ω β(x) ⋅ w(x) dx.
Since A q ′ is closed and convex in (W 1,q ′ (Ω)) * × (L q ′ (Ω)) d , assertion (i) follows. Using that B q (m, w) < +∞ implies that w(x) = 0 a.e. where m(x) = 0, the definition of A q ′ implies that (2.17)

sup (α,β)∈A q ′ Ω m(x) dα(x) + Ω β(x) ⋅ w(x) dx = sup β∈L q ′ (Ω) d Ω [-H(x, -β(x))m(x) + β(x) ⋅ w(x)] dx, = sup β∈L q ′ (Ω) d Ω [-H(x, -β(x)) + β(x) ⋅ v(x)] m(x) dx.
An element β ∈ L q ′ (Ω) d maximizes the above expression iff for a.e. x ∈ {m > 0} we have that

v(x) = -∇H(x, -β(x)). Since β ∈ L q ′ (Ω) d , if v ∉ L q (Ω) d Lemma 2
.1 implies that the previous relation cannot be satisfied and so ∂B q (m, w) = ∅. On the other hand, if v ∈ L q (Ω) d then any β ∈ L q ′ (Ω) d satisfying a.e. in {m > 0} that β(x) = -∇H * (x, -v(x)) optimizes the last expression in (2.17). Therefore, using the definition of A q ′ we readily get that if (α, β) ∈ ∂B q (m, w), then

β(x) = -∇H * (x, -v(x)) and α {m > 0} = -H(∇H * (x. -v)) a.e. in {m > 0}.
The result follows.

Remark 2.2. A generalization of the previous result to the case when 1 < q ≤ d could be interesting by extending the techniques in [START_REF] Brézis | Intégrales convexes dans les espaces de Sobolev[END_REF]. However, since our results in the next section are intrinsically related to the assumption q > d, we have preferred to provide a direct and self-contained proof in this case.

The variational problems

Let us fix q > d. In order to define the variational problems we are interested in, we introduce first the data and our assumptions. Let

W 1,q + (Ω) ∶= {m ∈ W 1,q (Ω) ; m(x) ≥ 0 ∀ x ∈ Ω},
and

F ∶ W 1,q (Ω) → R be such that (3.1)
F is weakly lower semicontinuous, Gâteaux-differentiable in W 1,q + (Ω) and satisfies that ∀ R > 0

; ∃ C R ∈ R such that if m ∈ W 1,q (Ω) and 0 ≤ m(x) ≤ R ∀ x ∈ Ω, then F(m) ≥ C R .
Given w ∈ L q (Ω) d let us consider the following elliptic PDE, with Neumann boundary conditions, (3.2)

-∆m + div(w) = 0, in Ω, (∇mw) ⋅ n = 0, on ∂Ω.

We say that m ∈ W 1,q (Ω) is a weak solution of (3.2) if

(3.3) Ω ∇m(x) ⋅ ∇ϕ(x) dx = Ω w(x) ⋅ ∇ϕ(x) dx ∀ φ ∈ C 1 (Ω).
Let us use the notation Y ∶= {ℓ ∈ W 1,q ′ (Ω) ; ∫ Ω ℓ(x) dx = 0}. Since, div(w) ∈ Y * , the results in [GM12, Section 7.1] and [START_REF] Mészáros | A variational approach to second order mean field games with density constraints: the stationary case[END_REF]Appendix] imply the existence of a unique m ∈ W 1,q (Ω) such that m is a weak solution of (3.3) and ∫ Ω m(x) dx = 1. Moreover, there exists C > 0 such that

(3.4) ∇m q ≤ C w q ,
and so, by the Poincaré-Wirtinger inequality, there exists C > 0 such that (3.5) m 1,q ≤ C ( w q + 1) .

Note that (3.3) can be written as Am

+ Bw = 0, with A ∶ W 1,q (Ω) → Y * and B ∶ L q (Ω) d → Y * being linear bounded operators defined as ⟨Am, φ⟩ Y * ,Y ∶= Ω ∇m ⋅ ∇φ dx, ⟨Bw, φ⟩ Y * ,Y ∶= - Ω w ⋅ ∇φ dx ∀ φ ∈ Y,
where ⟨⋅, ⋅⟩ Y * ,Y denotes the duality product between Y * and Y . Now, let us define

J ∶ W 1,q (Ω)× L q (Ω) d → R and G ∶ W 1,q (Ω) × L q (Ω) d → Y * × R as J (m, w) ∶= B q (m, w) + F(m), G ∶= (G 1 , G 2 ) ∶= Am + Bw, Ω m(x) dx -1 .
The first variational problem we consider is

(P 1 ) inf m∈W 1,q (Ω), w∈L q (Ω) d J (m, w) such that G(m, w) = 0.
In the second variational problem we impose a density constraint: let κ ∈ W 1,q (Ω) be such that

(3.6) κ ∶= min x∈Ω κ(x) > 0 and Ω κ(x) dx > 1.
Given a representative of κ, still denoted by κ, we define the set

C ∶= {m ∈ C(Ω) ; m(x) ≤ κ(x) ∀ x ∈ Ω}.
We consider the problem

(P 2 ) inf m∈W 1,q (Ω), w∈L q (Ω) d J (m, w) such that G(m, w) = 0, m ∈ C.
Note that q > d and the Sobolev embeddings imply that W 1,q (Ω) ↪ C(Ω) and so the constraint m ∈ C in (P 2 ) is well-defined.

Remark 3.1. Theorem 2.4(i) and (3.1) imply that the cost functional J in (P 1 ) and (P 2 ) is weakly lower semicontinuous. On the other hand, it is not necessarily convex.

3.1. Existence of solutions of the variational problems. In this subsection we prove that both problems (P 1 ) and (P 2 ) admit at least one solution. The proof of existence of solutions of problem (P 2 ) follows the same lines than the proof of [MS15, Theorem 3.1]. The proof of existence of solutions for problem (P 1 ) introduces an artificial density constraint and uses the existence of solutions of (P 2 ). Given a Lebesgue measurable set A ⊆ R d we denote by A its Lebesgue measure.

Theorem 3.1. Assume that (3.1) holds. Then, problem (P 2 ) has at least one solution (m, w). If in addition, F is bounded from below in W 1,q + by a constant C F ∈ R, then problem (P 1 ) also admits at least one solution (m, w). Moreover, in the latter case,

(3.7) m ∞ ≤ max 2c 0 c 1 , (2c 0 c 1 ) q qC q-1 1 (F(1 Ω ) + 2C 2 -C F ) , and
w q q ≤ qC q-1 1 (F(1 Ω ) + 2C 2 -C F ) m q-1 ∞ ,
where C 1 and C 2 satisfy (2.3) and c 0 and c 1 depend only on the geometry of Ω.

Proof. We first prove the assertion for problem (P 2 ), where the density constraint allows to obtain directly some bounds on any minimizing sequence. Define

(3.8) κ ∶= κ κ 1 and so 0 < κ(x) < κ(x) ∀ x ∈ Ω and Ω κ(x) dx = 1.
By the surjectivity result in [MS15, Lemma A.1], there exists ŵ ∈ L q (Ω) d such that Aκ + B ŵ = 0. By (2.4), we have that

B q (κ, ŵ) ≤ C q-1 1 q Ω ŵ(x) q κ q-1 dx + C 2 + F(κ) < +∞,
which implies that the infimum in (P 2 ) is not +∞. Now, let (m n , w n ) be a minimizing sequence and set κ ∶= max x∈Ω κ(x). The previous discussion implies the existence of C > 0 such that J (m n , w n ) ≤ C for all n ∈ N. In particular, by (2.2), w n = 0 a.e. on the set {m n = 0}. Since (m n , w n ) is feasible, we get that (1 m q-1 n )½ {mn >0} ≥ κ 1-q and so (2.4) and (3.1), with R = κ, imply that (3.9)

Ω w n q dx = {mn>0} w n q dx ≤ q(C 1 κ) q-1 (C + C 2 -C κ ) .
Therefore, the sequence w n is bounded in L q (Ω) d and so, by (3.5), the sequence m n is bounded in W 1,q (Ω). Therefore, extracting a subsequence, we obtain the existence of (m, w) such that m n converges weakly to m in W 1,q (Ω) and w n converges weakly to w in L q (Ω) d . By passing to the weak limit, we get that G(m, w) = 0 and m ∈ C. Finally, using that J is weakly lower semicontinuous we get that (m, w) solves (P 2 ). Now, let us prove existence for (P 1 ) under the additional assumption on the boundedness from below of F in W 1,q + (Ω). Let γ > 1 Ω . Since ( m, ŵ) ∶= (1 Ω , 0) is feasible for problem (P 2 ), with κ(x) ≡ γ, we have that (P 2 ) admits at least one solution. We will show that any such solution (m γ , w γ ) satisfies that m γ ∞ ≤ κ for some constant κ > 0 which is independent of γ. This will prove the result since any solution (m κ, w κ) of (P 2 ) with κ(x) ≡ κ solves (P 1 ). Indeed, if there is a feasible (m, w) for problem (P 1 ) such that J (m, w) < J (m κ, w κ) then since there exists κ ′ > 0 such that m ∞ ≤ κ ′ (because m ∈ W 1,q (Ω)) we have that J (m κ ′ , w κ ′ ) ≤ J (m, w), where (m κ ′ , w κ ′ ) is a solution of (P 2 ) with κ ≡ κ ′ , and m κ ′ ≤ κ which contradicts the optimality of (m κ, w κ) in (P 2 ) with κ = κ.

Let us denote by c 0 > 0 a constant such that m ∞ ≤ c 0 m 1,q for all m ∈ W 1,q (Ω) and by c 1 > 0 the constant in (3.5). Thus, any solution of G(m, w) = 0 satisfies that m ∞ ≤ c 0 c 1 ( w q + 1) ≤ 2c 0 c 1 max{ w q , 1}. Now, let us fix a solution (m γ , w γ ) of (P 2 ) with κ ≡ γ. If w γ q ≤ 1, then m γ ∞ ≤ 2c 0 c 1 , so let us assume that w γ q > 1 and so m γ ∞ ≤ 2c 0 c 1 w γ q . Since w γ vanishes a.e. in {m γ = 0}, arguing as in (3.9) we get that (3.10)

w γ q q ≤ qC q-1 1 (B q (1 Ω , 0) + F(1 Ω ) + C 2 -C F ) (2c 0 c 1 w γ q ) q-1 .
Since (2.4) implies that B q (1 Ω , 0) ≤ C 2 , we find that

w γ q ≤ qC q-1 1 (F(1 Ω ) + 2C 2 -C F ) (2c 0 c 1 ) q-1
and so m γ ∞ ≤ (2c 0 c 1 ) q qC q-1

1 (F(1 Ω ) + 2C 2 -C F ) .
The result follows.

3.2. Existence of solutions of the Mean Field Game systems. We recall that for a nonempty closed and convex set K, the normal cone N K (x) to K at x is defined as N K (x) ∶= ∂χ K (x).

We have the following existence result.

Theorem 3.2. Assume that (3.1) holds and that F is bounded from below in W 1,q + (Ω). Then, there exists (m, u, λ)

∈ W 1,q (Ω) × Y × R such that (3.11) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ A * u + H(⋅, ∇u) + λ = DF(m), -∆m -div m∇ ξ H(⋅, ∇u) = 0, in Ω, (∇m + m∇ ξ H(⋅, ∇u)) ⋅ n = 0, on ∂Ω, Ω m dx = 1, m(x) > 0, in Ω,
where the second equation, with its boundary condition, is satisfied in the weak sense (see (3.3)).

Proof. Theorem 3.1 yields the existence of a solution (m, w) ∈ W 1,q (Ω) × L q (Ω) d of problem (P 1 ). By definition,

Bq (m, w) + F(m) ≤ Bq ( m, w) + F( m) ∀ ( m, w) ∈ W 1,q (Ω) × L q (Ω) d ,
where Bq (m, w) ∶= B q (m, w) + χ G -1 (0) (m, w). Since Bq (m, w) is finite, the Gâteaux differentiability of F and the convexity of Bq imply that

-⟨DF(m), m⟩ ≤ lim inf τ ↓0 Bq (m + τ m, w + τ w) -Bq (m, w) τ = inf τ >0 Bq (m + τ m, w + τ w) -Bq (m, w) τ ,
where we have denoted by ⟨⋅, ⋅⟩ the duality product between (W 1,q (Ω)) * and W 1,q (Ω). Taking τ = 1 in the last term of the previous inequality, we get that

Bq (m, w) -⟨DF(m), m⟩ ≤ Bq (m + m, w + w) ∀ ( m, w) ∈ W 1,q (Ω) × L q (Ω) d , from which (-DF(m), 0) ∈ ∂ Bq (m, w). Set ( m, ŵ) ∶= (1 Ω , 0). Since χ G -1 (
0) ( m, ŵ) = 0 and B q is finite and continuous at ( m, ŵ) (since q > d), we have (see e.g. [ABM06, Theorem 9.5.4(b)])

(3.12) (-DF(m), 0) ∈ ∂ Bq (m, w) = ∂B q (m, w) + ∂χ G -1 (0) (m, w) = ∂B q (m, w) + N G -1 (0) (m, w).
In particular ∂B q (m, w) ≠ ∅ and so, by Theorem 2.4(ii), v = (w m)½ {m>0} ∈ L q (Ω) d . Lemma A.1, in the appendix, implies that m > 0 in Ω, hence Theorem 2.4(ii) implies that

∂B q (m, w) = {(-H(⋅, ∇H * (⋅, -v)), -∇H * (⋅, -v))} .
On the other hand, by [MS15, Lemma A.1] we have that G is surjective and so (see e.g. [START_REF] Bonnans | Perturbation analysis of optimization problems[END_REF])

N G -1 (0) (m, w) = {DG(m, w) * (û, λ) ; (û, λ) ∈ Y × R}, = {(A * û + λ, B * û) ; (û, λ) ∈ Y × R}.
Therefore, since B * û = -∇û, we get the existence of (û, λ) ∈ Y × R such that

(3.13) -A * û + H(⋅, ∇H * (⋅, -v)) -λ = DF(m), ∇H * (⋅, -v) = -∇û.
Thus, defining u = -û and λ = -λ we get the first equation in (3.11). On the other hand, since w = mv and v = -∇ ξ H(⋅, ∇u), the remaining equations in (3.11) also hold true.

The proof of the following result, concerning problem (P 2 ), is analogous to the previous one (see also Theorem 4.1 and Corollary 4.2 in [START_REF] Mészáros | A variational approach to second order mean field games with density constraints: the stationary case[END_REF]), hence we omit it. Notice, however, that the extra assumption on the global lower bound for F is not needed, since existence also holds true when we only assume (3.1) (see Theorem 3.1).

Theorem 3.3. Assume that (3.1) holds. Then, there exists (m, u, p, λ) ∈ W 1,q (Ω)×Y ×M (Ω)×R such that

(3.14) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ A * u + H(⋅, ∇u) -p + λ = DF(m), -∆m -div m∇ ξ H(⋅, ∇u) = 0, in Ω, (∇m + m∇ ξ H(⋅, ∇u)) ⋅ n = 0, on ∂Ω, Ω m dx = 1, 0 < m(x) ≤ κ(x), in Ω, spt(p) ⊆ {m = κ}, p ≥ 0, in Ω
where the second equation, with its boundary condition is satisfied in the weak sense (see (3.3)).

Remark 3.2. In the above theorem, the variable p plays the role of the Lagrange multiplier associated to the density constraint m ≤ κ.

We discuss now the uniqueness of solutions of systems (3.11) and (3.14) under a convexity assumption on F. Proposition 3.4. If F is strictly convex in W 1,q + (Ω) then the solutions of (3.11)-(3.14) are unique.

Proof. Let us consider first (3.11). Since F is convex, (P 1 ) is a convex problem and so if (m, u, λ) solves (3.11) then (m, -m∇ ξ H(⋅, ∇u)) solves (P 1 ). Thus, the uniqueness of m is a straightforward consequence of the strict convexity of F. Since m > 0 and for all x ∈ Ω the map

R d ∋ ξ ↦ ∇ ξ H(x, ξ) is injective and R d ∋ w ↦ H * x, w m(x) m(x)
is strictly convex we have that ∇u is unique. Thus, uniqueness of u in Y follows and, as a consequence, the first equation in (3.11) yields the uniqueness of λ. The proof of uniqueness of (m, u) for system (3.14) is the same as the previous one. By considering test functions supported in {0 < m < κ} we get the uniqueness of λ, from which the uniqueness of p follows. Let us detail the novelty of our results. Compared to [START_REF] Mészáros | A variational approach to second order mean field games with density constraints: the stationary case[END_REF], when q > d, we consider more general Hamiltonians and we prove that m is strictly positive in Ω which allows us to eliminate the Lagrange multiplier µ ∈ M (Ω) from the system (M F G) q in [START_REF] Mészáros | A variational approach to second order mean field games with density constraints: the stationary case[END_REF]. As we have seen, we can also get rid of the density constraint and prove, using variational methods, the existence of solutions of (3.11). Most importantly, we can consider, for both systems (3.11) and (3.14), rather general right-hand sides for the HJB equation, since we allow F to be non-convex. As an example of a class of functions we can deal with in (3.11), we consider

(3.15) F(m) ∶= Ω F (x, m(x), ∇m(x)) dx,
where

F ∶ Ω × R × R d → R is a Carathéodory function, i.e. for all (z, ξ) ∈ R × R d we have F (⋅, z, ξ)
is measurable and for a.e. x ∈ Ω the function F (x, ⋅, ⋅) is continuous. In addition, suppose that: (i) For a.e. x ∈ Ω and all z ∈ R the function F (x, z, ⋅) is convex.

(ii) For all R > 0 there exists γ ∈ L 1 (Ω) such that for a.e. x ∈ Ω (3.16)

F (x, z, ξ) ≥ γ(x) ∀ 0 ≤ z ≤ R, ξ ∈ R d .
(iii) For a.e. x ∈ Ω the function F (x, ⋅, ⋅) is differentiable. Moreover, for all R > 0 there exists

a 0 ∈ L 1 (Ω) and b 0 = b 0 (R) ≥ 0 such that for a.e. x ∈ Ω F (x, z, ξ) ≤ a 0 (x) + b 0 ξ q ∀ z ≤ R, ξ ∈ R d .
(iv) For all R > 0 there exists a 1 ∈ L 1 (Ω), a 2 ∈ L q ′ (Ω) and b 1 = b 1 (R) ≥ 0 such that for a.e.

x ∈ Ω, z ≤ R and ξ ∈ R d we have that

∂ z F (x, z, ξ) ≤ a 1 (x) + b 1 ξ q , ∇ ξ F (x, z, ξ) ≤ a 2 (x) + b 1 ξ q-1 .
Since q > d, assumptions (i)-(ii) imply the weak lower semi-continuity of F (see e.g. [Dac08, Corollary 3.24]), conditions (ii)-(iii) imply that F(m) ∈ R for all m ∈ W 1,q (Ω) and the last condition (iv) implies that F is Gâteaux differentiable at every m ∈ W 1,q (Ω) (see e.g. the proof of [Dac08, Theorem 3.37]). Thus, assumption (3.1) for F is satisfied in this case. The Gâteaux derivative of F at m ∈ W 1,q (Ω) is given by linear continuous functional

⟨DF(m), z⟩ = Ω ∂ z F (x, m(x), ∇m(x)) z(x) + ∇ ξ F (x, m(x), ∇m(x)) ⋅ ∇z(x) dx,
for all z ∈ W 1,q (Ω). We obtain the following corollary of Theorem 3.2 and of Theorem 3.3 when F is independent of ∇m.

Corollary 3.5. Let f ∶ Ω × R → R be a Carathéodory function. Assume that for all R > 0 there exist a, γ 1 ∈ L 1 (Ω) such that for a.e. x ∈ Ω (3.17)

∫ z 0 f (x, z ′ ) dz ′ ≥ γ 1 (x) and f (x, z) ≤ a(x) ∀ z ≤ R.
Then, (MFG 2 ) admits at least one solution (m, u, p, λ)

∈ W 1,q (Ω) × Y × M (Ω) × R. If, in addition, there exists γ 2 ∈ L 1 (Ω) such that ∫ z 0 f (x, z ′ ) dz ′ ≥ γ 2 (x)
for a.a. x ∈ Ω and z ≥ 0, then system (MFG 1 ) admits at least one solution (m, u, λ) ∈ W 1,q (Ω) × Y × R. In both cases, we have the additional regularity u ∈ W 1,s (Ω) for all s ∈ (1, d (d -1)).

Proof. Consider F (x, z) ∶= ∫ z 0 f (x, z ′ ) dz ′ in (3.15). Since (3.17) implies (3.1), existence of a solution (m, u, p, λ) ∈ W 1,q (Ω) × Y × M (Ω) × R of (MFG 2 ) follows directly from Theorem 3.3. Analogously, ∫ z 0 f (x, z ′ ) dz ′ ≥ γ 2 (x) for all m ≥ 0 implies that F has a global lower bound on W 1,q + (Ω). Hence, existence of a solution (m, u, λ) ∈ W 1,q (Ω) × Y × R of (MFG 1
) follows directly from Theorem 3.2. Since in both systems m ∈ C(Ω), assumption (3.17) also implies that f (x, m(x)) ∈ L 1 (Ω) and the W 1,s regularity for u (s ∈ ( 1, d (d -1))), in both systems, follows from [Sta65, Théorème 9.1]. Now, we comment on some other possible choices of F. Remark 3.4. (i) A simple and interesting example is given by f (x, z) = z r , where r > 0. In this case, F (x, z) = 1 r+1 z r+1 if z > 0 and F (x, z) = -1 r+1 z r+1 , otherwise. Thus, the assumptions of Theorem 3.5 and Proposition 3.4 are satisfied and so systems (MFG 1 ) and (MFG 2 ) admits unique solutions. Note that the growth of f is arbitrary, showing one the advantages of the variational approach (compare with [GPV14, [START_REF] Cirant | Multi-population mean field games systems with Neumann boundary conditions[END_REF][START_REF] Pimentel | Regularity for second order stationary mean-field games[END_REF] where the growth of f is restricted). On the other hand, if we consider f (x, z) =z r , its primitive is not bounded from below in [0, ∞) and the assumptions of Theorem 3.5 are not satisfied for problem (MFG 1 ). However, they are satisfied for (MFG 2 ) and the existence of at least one solution to (MFG 2 ) is ensured also in this case.

(ii) In order to exemplify the possible dependence of F on ∇m, let us take F(m) ∶= 1 2 ∫ Ω ∇m 2 dx. In this case Theorem 3.2 yields the existence of weak solution of

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ -∆u + H(⋅, ∇u) + λ = -∆m, in Ω, ∇(u -m) ⋅ n = 0, on ∂Ω, -∆m -div m∇ ξ H(⋅, ∇u) = 0, in Ω, (∇m + m∇ ξ H(⋅, ∇u)) ⋅ n = 0, on ∂Ω, Ω m dx = 1, m(x) > 0, in Ω.
Moreover, by Proposition 3.4, the solution is unique.

(iii) We can also consider a non-local dependence on m. For instance, setting Ω ε ∶= {x ∈ Ω ; d(x, ∂Ω) > ε}, where d(⋅, ∂Ω) is the distance function to ∂Ω, we can take F as before and

F(m) ∶= Ωε F (x, ρ 0 * m, ρ 1 * ∂ x 1 m(x), ..., ρ d * ∂ x d m(x)) dx,
for some given regular kernels ρ 0 , ..., ρ d supported on B ε (0). In this case, it is easy to check that (3.1) holds without requiring the convexity of F (x, z, ⋅).

Now let f ∶ Ω × R → R be a Carathéodory function satisfying that there exists γ 2 ∈ L 1 (Ω) such that ∫ z 0 f (x, z ′ ) dz ′ ≥ γ 2 (x) for a.e.
x ∈ Ω and z ≥ 0. Moreover, assume that for all R > 0 there exists a ∈ L 1 (Ω) such that for a.e. x ∈ Ω and z ≤ R we have that f (x, z) ≤ a(x). Under these assumptions, Corollary 3.5 ensures the existence of at least one solution (m, u, λ) ∈ W 1,q (Ω) × Y × R of (MFG 1 ). Using a bootstrapping argument, we show in the next result some additional local regularity properties for (m, u). Proposition 3.6. Consider the above setting and suppose that Ω ∋ x ↦ f (x, m(x)) ∈ R belongs to L r (Ω) for some r > d and that ∇ ξ H(x, ⋅) is Hölder continuous, uniformly on x ∈ Ω. Then, there exist α 0 , α 1 ∈ (0, 1) such that

(3.18) u ∈ C 1,α 0 loc (Ω) and m ∈ C 1,α 1 loc (Ω). Proof.
Step 1. We show that there exists k > d such that u ∈ W 2,k loc (Ω). By the classical Sobolev embeddings, this implies that u ∈ C 1,α 0 loc (Ω) (for some α 0 ∈ (0, 1)). Let r 1 ∈ (q ′ , d (d -1)). By Corollary 3.5, we have that ∇u q ′ ∈ L r 1 q ′ (Ω) and so, by (2.3), we have that H(⋅, ∇u) ∈ L r 1 q ′ (Ω).

Furthermore, since f (⋅, m(⋅)) ∈ L r (Ω) and r > δ 1 ∶= r 1 q ′ , the classical regularity theory for elliptic equations (see [START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF]) implies that u ∈ W 2,δ 1 loc (Ω). In particular, the Sobolev inequality (see e.g. [START_REF] Adams | Sobolev spaces[END_REF]) yields u ∈ W 1,

dδ 1 d-δ 1 loc
(Ω) and so ∇u q ′ ∈ L δ 2 loc (Ω) with δ 2 ∶= dδ 1 q ′ (d-δ 1 ) . We easily check that δ 2 > δ 1 and so we improve the regularity of u to obtain that u ∈ W 2,min{r,δ 2 } loc (Ω).

Thus, if δ 2 > d we obtain the first relation in (3.18). Otherwise, for i ≥ 2, inductively we define the sequence δ i+1 ∶= dδ i (d-δ i )q ′ . Since δ i+1δ i ≥ (q ′ + ddq ′ ) (dδ i )q ′ and q ′ + ddq ′ > 0, we get that δ i+1δ i ≥ (q ′ + ddq ′ ) dq ′ if d > δ i . Therefore, after a finite number of steps we get the existence of i * ≥ 2 such that δ i * > d and u ∈ W 2,k loc (Ω) with k ∶= min{r, δ i * } > d.

Step 2. Let us prove that m ∈ C 1,α 1 loc (Ω) for some α 1 ∈ (0, 1). Since m ∈ W 1,q (Ω) and q > d, we already have that m is Hölder continuous. Having u ∈ C 1,α 0 loc (Ω), this implies that ∇u ∈ C 0,α 0 loc (Ω) d , hence m∇ ξ H(⋅, ∇u) ∈ C 0, α loc (Ω) d , for some α ∈ (0, 1). Using a Schauder-type estimate (see [START_REF] Giaquinta | An introduction to the regularity theory for elliptic systems, harmonic maps and minimal graphs[END_REF]Theorem 5.19]) we get that m ∈ C 1,α 1 loc (Ω) for some α 1 ∈ (0, 1). Remark 3.5. If f satisfies (3.17), with a ∈ L r (Ω) for some r > d, and (m, u, p, λ) solves (MFG 2 ), we have that (u, m) admits the regularity (3.18) locally in the open set {0 < m < κ} (see [START_REF] Mészáros | A variational approach to second order mean field games with density constraints: the stationary case[END_REF]Proposition 4.3] for a similar result in a simpler case).

An application to multipopulation systems

In this section we show a simple application of our results to the study of systems of the form

(MFG N ) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ -∆u i + H i (⋅, ∇u i ) + λ i = f i (x, (m i ) N i=1 ), in Ω, ∇u i ⋅ n = 0, on ∂Ω, -∆m i -div m i ∇ ξ H i (⋅, ∇u i ) = 0, in Ω, (∇m i + m i ∇ ξ H i (⋅, ∇u i )) ⋅ n = 0, on ∂Ω, Ω m i dx = 1, m i (x) > 0, in Ω, i = 1, ..., N,
where N ≥ 2 (see [START_REF] Cirant | Multi-population mean field games systems with Neumann boundary conditions[END_REF][START_REF] Bardi | Nonlinear elliptic systems and mean-field games[END_REF]). Here (m i ) N i=1 describes the densities of N populations. The Hamiltonians H i ∶ Ω × R d → R (i = 1, . . . , N ) are supposed to satisfy the assumptions in Section 2 (see in particular (2.3)). The given functions

f i ∶ Ω × [0, +∞) N → R (i = 1, . . . , N ) are such that for all ζ ∈ R N the function f i (⋅, ζ) is measurable and for a.e. x ∈ Ω the function f i (x, ⋅) is continuous. Suppose that (4.1) ∃ γ i ∈ L 1 (Ω) such that ∫ z 0 f i (x, z i , (ζ j ) j≠i ) dz i ≥ γ i (x) for a.e. x ∈ Ω, ∀ z ≥ 0, ∀ (ζ 1 , ..., ζ i-i , ζ i+1 , ..., ζ N ) ∈ [0, +∞) N -1 ,
where we have denoted

f i (x, z i , (ζ j ) j≠i ) ∶= f i (x, ζ 1 , ..., ζ i-i , z i , ζ i+1 , ..., ζ N ).
Moreover, we assume that for all i = 1, . . . , N

(4.2) ∀ R > 0, ∃ a i ∈ L 1 (Ω) such that f i (x, z, (ζ j ) j≠i ) ≤ a i (x), for a.e. x ∈ Ω, ∀ z ≤ R, ∀ (ζ 1 , ..., ζ i-i , ζ i+1 , ..., ζ N ) ∈ [0, +∞) N -1 ,
and that

(4.3) ∀ (ζ 1 , ..., ζ i-i , ζ i+1 , ..., ζ N ) ∈ [0, +∞) N -1 the map z ∈ [0, +∞) → f i (x, z, (ζ j ) j≠i ) ∈ R is non-decreasing.
We have the following result:

Proposition 4.1. Suppose that for all i = 1, . . . , N the function f i satisfies (4.1), (4.2) and (4.3). Then, system (MFG N ) admits at least one solution m = (m 1 , ..., m N ), u = (u 1 , ..., u N ) and λ = (λ 1 , ..., λ N ), where, for all i = 1, ..., N , m i ∈ W 1,q (Ω) and u

i ∈ W 1,s (Ω) (for all s ∈ (1, d (d -1)). Proof. Let us define F i ∶ W 1,q (Ω) N → R as F i (m i , (m j ) j≠i ) ∶= ∫ Ω F i (x, m i (x), (m j (x)) j≠i ) dx where F i (x, z, (ζ j ) j≠i ) ∶= ∫ z 0 f i (x, z i , (ζ j ) j≠i ) dz i .
Consider the space X ∶= W 1,q (Ω) × L q (Ω) d N , endowed with the weak-topology, and the setvalued map T ∶ X → 2 X defined as

T (x) ∶= N i=1 argmin (m ′ ,w ′ )∈W 1,q (Ω)×L q (Ω) d B q (m ′ , w ′ ) + χ G -1 (0) (m ′ , w ′ ) + F i (m ′ , (m j ) j≠i ) ,
where x = ((m 1 , w 1 ), ..., (m N , w N )). By Theorem 3.1, the embedding W 1,q (Ω) ↪ C(Ω) and our assumptions, we have that T (x) is a non-empty compact set for all x ∈ X. Assumption (4.3) implies that T (x) is also convex. Moreover, by (3.7) and (4.2), for x ∈ X we have the existence of c > 0, independent of x, such that mi 1,q ≤ c and wi q ≤ c for all m = ( m1 , ..., mN ) and w = ( w1 , ..., wN ) such that ( m, w) ∈ T (x). Therefore, defining

K c ∶= {(m, w) ∈ X ; m i 1,q ≤ c, w i q ≤ c, ∀ i = 1, ..., N }, we have that T (K c ) ⊆ K c . Now, let us prove that T is upper-semicontinuous, i.e. T -1 (M ) ∶= {x ′ ∈ X ; T (x ′ )∩M ≠ ∅} is closed for all closed sets M ⊆ X. Indeed, let x n = ((m n 1 , w n 1 ), ...., (m n N , w n N )) ∈ T -1 (M ) such that x n → x = ((m 1 , w 1 ), ...., (m N , w N ))
. By definition, we have the existence of xn = (( mn 1 , wn 1 ), ...., ( mn

N , wn N )) ∈ M such that B q ( mn i , wn i ) + χ G -1 (0) ( mn i , wn i ) + F i ( mn i , (m n j ) j≠i ) ≤ B q (m ′ , w ′ ) + χ G -1 (0) (m ′ , w ′ ) + F i (m ′ , (m n j ) j≠i ) ∀ (m ′ , w ′ ) ∈ W 1,q (Ω) × L q (Ω) d .
By (3.7) we have that xn is bounded in W 1,q (Ω) × L q (Ω) d N and so, up to some subsequence, there exists x = (( m1 , w1 ), ...., ( mN , wN )) such that xn → x in X and so, since M is closed, x ∈ M . Under our assumptions, the Lebesgue's dominated convergence theorem implies the weak continuity of F i in W 1,q (Ω) N , and so we can pass to the limit to obtain that x ∈ T -1 (M ). By Kakutani fixed-point theorem, there exists x = ((m 1 , w 1 ), ..., (m N , w N )) such that x ∈ T (x).

The result follows from Corollary 3.5.

Remark 4.1. (i) As we pointed out, the result in Proposition 4.1 is a simple consequence of the variational method we presented in the previous sections. We refer the reader to [Cir15, CV16, BF16, ABC17] for a more detailed study, and sharper results, based on PDE arguments tackling directly system (MFG N ).

(ii) The local regularity results presented in Proposition 3.6 for the one-population case directly extend to the solutions of system (MFG N ).

We can also consider the instance of (MFG N ) where the functions f i (i = 1, ..., N ) satisfy that there exists a Carathéodory function function F ∶ Ω × R N → R such that for a.e. x ∈ Ω the function F (x, ⋅) is differentiable and for all i = 1, ..., N we have that As suggested in [START_REF] Cirant | Multi-population mean field games systems with Neumann boundary conditions[END_REF]Remark 15], in this case system (MFG N ) can be found as the optimality condition of the optimization problem (P N ) inf m i ∈W 1,q (Ω), w i ∈L q (Ω) d i=1,...,N

N i=1

B q (m i , w i ) + χ G -1 (0) (m i , w i ) + Ω F (x, m 1 (x), ..., m N (x)) dx.

Indeed, suppose that F satisfies that there exists γ ∈ L 1 (Ω) such that Moreover, suppose that for all R > 0 there exists a ∈ L 1 (Ω) such that where we have used the fact that both ϕ and ϕ (restricted to U ) are admissible test functions in (3.3), a change of variable in the integrals and the properties that we have shown for S.

The regularity of v and [BKRS15, Corollary 1.7.2] imply that m(x) > 0. The result follows.

  for a.e. x ∈ Ω , and denote by M + (Ω) and M -(Ω) the subsets of nonnegative and nonpositive finite Radon measures of M (Ω), respectively. For a set D, we denote by ½ D its indicator function, i.e. ½ D (y) = 1 if y ∈ D and ½ D (y) = 0 otherwise.

Remark 3. 3 .

 3 The previous uniqueness result for solutions of (3.14) improves the one stated in [MS15, Remark 4.2].

  (4.4)f i (x, ζ i , (ζ j ) j≠i ) = ∂ ζ i F (x, ζ) for a.e. x ∈ Ω, ∀ ζ = (ζ 1 , ..., ζ N ) ∈ R N .

F

  (x, ζ) ≥ γ(x) for a.a. x ∈ Ω and for all ζ ∈ R N .

F

  (x, ζ) ≤ a(x) for a.a. x ∈ Ω and ζ ∈ R N such that ζ i ≤ R for all i = 1, ..., N .Then, arguing as in the proof of Theorem 3.1, we get the existence of a solution m = (m 1 , ..., m N ), w = (w 1 , ..., w N ) of (P N ), and so, mimicking the proof of Theorem 3.2, we get the following result:Defining T ∶ B d-1r (0) × (-r, r) → G as T (y, s) ∶= (y, ψ(y) + s), we have that T is a bi-Lipschitz map (an injective Lipschitz continuous map whose inverse is also Lipschitz continuous) andDT (y, s) = I d-1 0 ∇ψ(y) ⊺ 1 and DT (y, s) -1 = I d-1 0 -∇ψ(y) ⊺ 1 , for a.e. (y, s) ∈ B d-1 r (0) × (-r, r) (I d-1 ∈ R (d-1)×(d-1) denotes the identity matrix). Define thereflection map S ∶ G → G as (A.1) S(T (y, s)) = T (y, -s),which clearly leaves the points on ∂Ω ∩ G invariant and satisfies that S(S(x)) = x for all x ∈ G.Differentiating both sides of (A.1) yields for a.e. (y, s)∈ B d-1 r (0) × (-r, r) DS(T (y, s))DT (y, s) = DT (y, -s) I d-1 0 0 ⊺ R d-1 -1, from where we getDS(T (y, s)) = DT (y, -s) I d-1 0 0 ⊺ R d-1 -1 DT (y, s) -1 = I d-1 0 2∇ψ(y) ⊺ -1 .Thus, det(DS(x)) = -1, DS(x) -1 = DS(x), for every x ∈ G, and DS(⋅) is bounded on G. Note that since DS(T (y, s)) does not depend on s we have that DS(Sx) = DS(x).Let us introduce some notations. We setU ∶= G ∩ Ω, V ∶= G ∖ Ω = S(U ) and for b ∶ U → R, we define b ∶ V → U as b(x) ∶= b(S(x)) and b(x) = b(x), x ∈ U, b(x), x ∈ V, the a.e. extension of b to G (the definition on ∂Ω∩G is irrelevant since this set is L d -negligible). Arguing as in [Nit11, Lemma 3.3], if b ∈ W 1,q (U ), we have that b ∈ W 1,q (G), ∇ b = ∇b½ U + ∇b½ V and ∇b(x) = DS(x)∇b(S(x)) for a.e. x ∈ G. Now, since (m, v) satisfies (3.3) (with w = vm) for tests functions ϕ ∈ C 1 (Ω), defining v ∈ L q (G) d as v(x) = v(x), x ∈ U, DS(x)v(x), x ∈ V,(where v is understood componentwise), the pair ( m, v) solves a similar equation on G with test functionsϕ ∈ C 1 c (G). Indeed, let us take ϕ ∈ C 1 c (G) and compute (A.2) ∫ G (∇ m ⋅ ∇ϕmv ⋅ ∇ϕ) dx = ∫ U (∇ m ⋅ ∇ϕmv ⋅ ∇ϕ) dx + ∫ V (∇ m ⋅ ∇ϕmv ⋅ ∇ϕ) dx = ∫ U (∇m ⋅ ∇ϕmv ⋅ ∇ϕ) dx + ∫ V (∇m(x) ⋅ ∇ϕ(x)m(x)∇ϕ(x) ⊺ DS(x)v(x)) dx = ∫ V ∇ϕ(x) ⊺ DS(x)∇m(S(x)) dx -∫ V m(S(x))∇ϕ(x) ⊺ DS(x)v(S(x)) dx = ∫ U ∇ϕ(S(y)) ⊺ DS(S(y))∇m(y) dy -∫ U m(y)∇ϕ(S(y)) ⊺ DS(S(y))v(y) dy = ∫ U (∇m(y)m(y)v(y)) ⋅ ∇ϕ(y) dy = 0,

Acknowledgements

The second author was partially supported by the ANR (Agence Nationale de la Recherche) project ANR-16-CE40-0015-01 and by the Gaspar Monge Program for Optimization and Operation Research (PGMO) via the project PASTOR. Both authors were partially supported by the PGMO via the project VarMFGPDE. The authors also thank M. Cirant for suggesting them to look carefully whether some Harnack inequality holds in their framework. This observation led the authors to find the useful references in [BKRS15] and improve the results in the present paper.

Proposition 4.2. Suppose that f i , i = 1, ..., N , satisfy (4.4), with F satisfying (4.5)-(4.5). Then, system (MFG N ) admits at least one solution m = (m 1 , ..., m N ), u = (u 1 , ..., u N ) and λ = (λ 1 , ..., λ N ), where, for all i = 1, ..., N , m i ∈ W 1,q (Ω) and u i ∈ W 1,s (Ω) (for all s ∈ (1, d (d -1)).

Note that (4.4) is restrictive. On the other hand, the previous result does not require the strong boundedness condition (4.2) and the monotonicity assumption (4.3). Moreover, this framework allows us to introduce density constraints of the form m ∈ K, where

We suppose that κ ∈ W 1,q (Ω) satisfies κ(x) > 0, for all x ∈ Ω, and the weights (α i ) N i=1 satisfy (4.7)

Condition (4.7) implies that if for all i = 1, ..., N we define mi ∶= κ κ 1 and ŵi ∶= ŵ, where ŵ is such that B ŵ = -Aκ κ 1 (we know that such ŵ exists by [MS15, Lemma A.1]), then m ∶= ( m1 , ..., mN ), ŵ ∶= ( ŵ1 , ..., ŵN ) are feasible for problem

and m is an interior point to the constraint m ∈ K, i.e. ∑ N i=1 α i mi (x) < κ(x) for all x ∈ Ω. Therefore, we can argue as in the proof of Theorem 3.1 to show the existence of at least one solution of (P ′ N ) and then, following the proof of Theorem 3.1 (see also the proofs of Theorem 4.1 and Corollary 4.2 in [MS15]), we get the following result: Proposition 4.3. Suppose that for all i = 1, ⋯, N the function f i satisfies the assumptions of Proposition 4.2. Moreover, assume that (4.7) holds. Then system

admits at least one solution m = (m 1 , ..., m N ), u = (u 1 , ..., u N ), λ = (λ 1 , ..., λ N ) and p, where, for all i = 1, ..., N ,

We prove in this appendix a lemma which is crucial to show the strict positivity of the densities in (MFG 1 ) and (MFG 2 ).

Lemma A.1. Assume that Ω ⊂ R d satisfies the assumptions at the beginning of Section 2, q > d and let v ∈ L q (Ω) d . Suppose that m ∈ W 1,q (Ω) is a weak solution of

Then, m(x) > 0 for all x ∈ Ω Proof. Since v ∈ L q (Ω) d (q > d), by the Harnack's inequality proved in [START_REF] Trudinger | Linear elliptic operators with measurable coefficients[END_REF] (see also [START_REF] Bogachev | Fokker-Planck-Kolmogorov equations[END_REF] Corollary 1.7.2]) we have that m(x) > 0 for all x ∈ Ω. It remains to study the positivity of m on the boundary ∂Ω. To achieve this, we use a standard reflection argument following [START_REF] Nittka | Regularity of solutions of linear second order elliptic and parabolic boundary value problems on Lipschitz domains[END_REF]. Let x ∈ ∂Ω. Since Ω is supposed to be Lipschitz domain, in a neighborhood of x the boundary ∂Ω can be represented as the graph of a Lipschitz function ψ ∶ R d-1 → R. Without loss of generality for these local considerations we may suppose that x = 0 and that

0)}, for some r > 0 small (otherwise one can perform some isometric transformations on Ω under which all the hypotheses that we assumed on the data remain invariant). We have denoted by B d r (z) and B d-1 r (z) the balls of radius r centered at z in R d and R d-1 , respectively. Consider a strip-like domain G ∶= (y, ψ(y) + s) ∶ y ∈ B d-1 r (0), s ∈ (-r, r) , such that Ω ∩ G = (y, ψ(y) + s) ∶ y ∈ B d-1 r (0), s ∈ (0, r) .