Denis Ullmo

Igor Swiecicki

Thierry Gobron

M Lions

R P Huang

P E Malhamé

Caines

Quadratic mean field games

Mean field games were introduced independently by J-M. Lasry and P-L

I. INTRODUCTION

Differential Games represent a field of mathematics at the frontier between optimization problems and Game Theory. It relates to optimization problems in the sense that they model socio-economic phenomena in which agents have control of some parameters (e.g. their velocity, or the amount of resource they dedicate to some goal), which can be modified in the course of time in order to minimize some given cost functional. It relates also to Game Theory because the cost function of a given agent is assumed to depend not only on his (or her) own state but also on the one of the other agents involved, implying a strategic approach to the parameter choices.

When the number of agents becomes large, differential games may become technically very complex. Practical applications for a large number of players are therefore extremely difficult to implement [START_REF] Khan | Non-cooperative games with many players[END_REF]. In some aspects, this situation is reminiscent of the problems encountered by physicists when considering systems with a large number of interacting components, for which an exact treatment is in most of the cases intractable, but where a "mean field" approximation provides in many circumstances a very decent description. Owing explicitely to this source of inspiration, P-L. Lions and J-M. Lasry [START_REF] Lasry | Jeux à champ moyen. II -Horizon fini et contrôle optimal[END_REF][START_REF] Lasry | Jeux à champ moyen. I -Le cas stationnaire[END_REF][START_REF] Lasry | Mean field games[END_REF], and independently M.Huang, R.P. Malhamé and P. E. Caines [START_REF] Huang | Large population stochastic dynamic games: closed-loop McKean-Vlasov systems and the Nash certainty equivalence principle[END_REF], introduced a "mean field" approximation scheme adapted to these differential games. The key idea at the root of Mean Field Games (MFG in the following), is that the very complexity associated with a large number of players allows for a drastic simplification when considering that a given agent is not really sensitive to the individual choices of every agent, but only to an averaged quantity (the mean field) which aggregates the decisions made by the other participants to the game.

Since the first articles in 2005-2006, Mean Field Games has developped into a very active field of research and undergone a rapid growth in several directions : On the more formal side, important results have been obtained on the existence and uniqueness of the associated PDEs ( [START_REF] Cardaliaguet | Notes on mean field games (from P.-L. Lions' lectures at Collège de France)[END_REF][START_REF] Gomes | Mean field games models -a brief survey[END_REF]) or on the convergence of a many player game to its mean field version ([9, 19, 22]). Significant work has been also done in order to find numerical schemes able to handle this kind of problems ([3, 36, 50]). More recently, application oriented studies have been performed in the context of finance ( [START_REF] Cardaliaguet | Mean field game of controls and an application to trade crowding[END_REF][START_REF] Lachapelle | Efficiency of the price formation process in presence of high frequency participants: a mean field game analysis[END_REF][START_REF] Lasry | Mean field games[END_REF]), economic problems ( [START_REF] Achdou | Partial differential equation models in macroeconomics[END_REF][START_REF] Achdou | A long-term mathematical model for mining industries[END_REF][START_REF] Guéant | Mean Field Games and Applications[END_REF]), or engineering [START_REF] Kizilkale | Chapter twenty -collective target tracking mean field control for markovian jump-driven models of electric water heating loads[END_REF][START_REF] Mériauxi | Mean field energy games in wireless networks[END_REF]. This recent period is thus marked by important results coming from the mathematical and engineering science communities, with more recent contributions from economists, but very little involvement from physicists.

One important goal of this review is to demonstrate that this has not to be the case, and that, as they did for other problems that have emmerged in the context of social sciences (cf eg [7, 14-16, 28, 29, 31]), physicists can bring an interesting and important point of view to the study of Mean Field Games. Indeed, although the MFG equations allow for a drastic simplification of the original problem, they are still very difficult to analyze. Very few exact solutions exist, mainly in the static case or in very simplified settings [START_REF] Almulla | Two numerical approaches to stationary mean-field games[END_REF][START_REF] Bardi | Explicit solutions of some linear-quadratic mean field games[END_REF][START_REF] Guéant | A reference case for mean field games models[END_REF][START_REF] Swiecicki | Phase diagram" of a mean field game[END_REF][START_REF] Laguzet | Individual vaccination as Nash equilibrium in a SIR model with application to the 2009-2010 influenza a (H1N1) epidemic in france[END_REF]. Furthermore the numerical schemes that have been developed, in spite of their quantitative accuracy, do not necessarily provide a complete understanding of the mechanism at work in Mean Field Games. A physicist's approach to MFG, in which ones develops concurrently both an intuition of the qualitative behaviors and quantitatively accurate approximation schemes, is therefore very much needed. Thus, rather than an exhausting review of existing mathematical results, the aim of this paper is to provide a good entry point for physicists to this fascinating new and rapidly developing field of Mean Field Games. For a recent survey of mathematical results, the reader is referred to the existing reviews [START_REF] Cardaliaguet | Notes on mean field games (from P.-L. Lions' lectures at Collège de France)[END_REF][START_REF] Gomes | Mean field games models -a brief survey[END_REF][START_REF] Guéant | Mean Field Games and Applications[END_REF] and books [START_REF] Bensoussan | Mean Field Games and Mean Field Type Control Theory[END_REF][START_REF] Carmona | Probabilistic Theory of Mean Field Games with Applications[END_REF][START_REF] Carmona | Probabilistic Theory of Mean Field Games with Applications[END_REF] covering in very much details the results obtained by the mathematical community.

In this paper, we shall on the other hand focus on a particular class of Mean Field Games, the so-called "quadratic Mean Field Games" which is still representative of the subject at large, and has the peculiarity of exhibiting a direct, formal, and deep relationship with the non-linear Schrödinger equation. This connection brings a wide set of existing tools to this field, which because of the long history of the non-linear Schrödinger equation in physics are both extremely familiar to physicists and lead to very significant progresses in the comprehension of the solutions of Mean Field Games models.

The body of this article is divided into five parts. In section II, we introduce in details the generic form of Mean Field Games, and present some of its extensions and applications. This section will be the only one that genuinely qualify as a [bird's-eye] survey of the mathematical literature.

The rest of the article is devoted to the study of the behavior of the quadratic mean field games, stressing in particular how tools which are well known to the physics community can be used effectively to understand the corresponding models.

In section III, we derive the connection between quadratic MFG systems and the non linear Schrödinger equation. As a first application, we present a few exact relations satisfied by some of the MFG statistics and revisit the stationary problem. We also derive exact solutions of "soliton"-type for some specific MFG problems, part of them being previously known, and some others new. In the two following sections, we study thoroughly the case of strong interactions, both in stationary and non-stationary settings. Section IV presents the main methods and results; in section V we extend these results in various ways (higher dimensions, general initial conditions, etc..).

Finally, in section VI we consider the case of weak interactions that can be cast into the same framework.

II. MEAN FIELD GAMES

In this section, we first describe in some details the construction of Mean Field Games in the form in which they were originally introduced ( [START_REF] Lasry | Jeux à champ moyen. II -Horizon fini et contrôle optimal[END_REF][START_REF] Lasry | Jeux à champ moyen. I -Le cas stationnaire[END_REF]), and review some generalizations, as well as some example of applications.

A. The "Mean Field Games" Paradigm

In its original form, a MFG problem can be set as follow. As a starting point, we consider a differential game with N players (or agents). Their individual states are described by a continuous variable X i ∈ R d , i = 1 . . . N , which, depending on the context, may represent a physical position, the amount of resources owned by a company, the house temperature in a network of controlled heaters, etc.. These state variables evolve according to some controlled dynamics, which here is assumed to be described a linear Langevin equation,

dX i t = a i t dt + σdW i t (2.1)
with initial conditions X i t = x i 0 . For simplicity, σ is chosen here to be a constant, each of the d components of W i is an independent white noise of variance one, and the control parameter is the velocity a i t . This control is adjusted in time by the agent i in order to minimize a cost functional over a time interval [0, T ]

c[a i ](x i t , t) = T t L(X i τ , a i τ ) -Ṽ [m τ ](X i τ ) dτ noise + c T (X i T ) noise . (2.2)
The cost functional is an average over all realizations of the noise for an agent starting at x i t at time t and represents the rational expectations of each player. In differential games, it is made of two parts: a time integral over a "running cost", here L(x, a) -Ṽ [m](x), and a "final cost", c T (x), depending on the state of the agent at the end of the optimization period T . The running cost here splits into two terms, a "free Lagrangian part", L(x, a), depending only on the agent's state x and control a, and a "potential" term, -V [m t ](x), which takes into account the interaction between agents and is a functional of the empirical density of agents m t in the state space,

m t (x) = 1 N j δ(x -X j (t)) . ( 2.3)
Implicitly, we have assumed that all agents have an identical behavior in the sense that they may differ only by their initial conditions and the subsequent choices of control parameters. Finally, as usual in control theory [START_REF] Bertsekas | Dynamic Programming and Optimal Control[END_REF] one introduces a value function which is defined as the minimum, over all controls, of the cost function, given the initial condition x at time t. It is independent on the agent label i and reads u t (x) ≡ min a c[a](x, t) .

(2.4)

Up to now, we have only described a continuous differential game, and it has to be made clear that even in the above simplified setting, there is no hope to solve this system beyond a very few number of players. The aim of Mean Field Games is to provide a frame in which such a system can be analyzed in the limit of a large number of players. One thus renounces to follow each agent individually, but rather describes the system at a statistical level through the density of agents m t (x). Now, the essential assumption of Mean Field Games is that the density of agents becomes deterministic in the large N limit. As a consequence, like in standard mean field theories, for a given time dependent density m t , the trajectories of the different players decouple, and the problem reduces to an optimization problem for a single agent. In particular, using Ito's calculus, the value function for each player can be shown to be solution of the following Hamilton-Jacobi-Bellman equation (we give a brief sketch of the derivation in appendix A)

       ∂ t u t (x) + H(x, ∇ ∇ ∇u t (x)) + σ 2 2 ∆u t (x) = Ṽ [m t ](x) u T (x) = c T (x) (2.5)
where

H(x, p) ≡ inf α α α (L(x, α α α) + p.α α α)).
The optimal control is then the value of a which realizes the above infimum, namely a = a * t (x) = ∂H ∂p (x, ∇u t (x)). Substituting into the Langevin equation (2.1), the optimized agent density m t then evolves from an initial condition m 0 according to the Kolmogorov equation [START_REF] Risken | The Fokker-Planck Equation: Methods of Solution and Applications[END_REF]

       ∂ t m t (x) + ∇ ∇ ∇.(m t (x) a * t (x)) - σ 2 2 ∆m t (x) = 0 m 0 (x) = m 0 (x) . ( 2.6)
By consistency, the two agent densities (i.e. the one used in the optimization leading to Eq. (2.5) and the one resulting from Eq. (2.6)) need to coincide and we are left with the set of partial differential equations (2.5) and (2.6) which are coupled together through the terms Ṽ [m t ] and a * t . They are both of diffusion type, but respectively backward and forward in time. Together, they form the system of equations defining the Mean field Game which is expected to describe the large N limit of our initial differential game [START_REF] Lasry | Jeux à champ moyen. II -Horizon fini et contrôle optimal[END_REF].

In the limit of a very large optimization period T → ∞, this system of equation has a remarkable property, proven under some specific conditions by P. Cardialaguet et al. [START_REF] Cardaliaguet | Long time average of mean field games with a nonlocal coupling[END_REF]. In a wide time span, when sufficiently far from both limits, t = 0 and t = T , the system stays in a permanent regime where the solution (m t (x), u t (x)) of the MFG system remains well approximated by the (m e (x), u e (x)e λet ), where the couple (m e (x), u e (x)) is solution of an analog ergodic problem (we assume here its existence and unicity)

         -λ e + H(x, ∇ ∇ ∇u e (x)) + σ 2 2 ∆u e (x) = Ṽ [m e ](x) ∇ ∇ ∇(m e (x)a e (x)) - σ 2 2 ∆m e (x) = 0 , ( 2.7)
where

H(x, p) = inf α (L(x, α) + p • α)), a e (x) = ∂H ∂p (x, ∇ ∇ ∇u e (x)
) and λ e an appropriately chosen constant. This system of equations is referred to as the stationary Mean Field Games system of equations. In the sequels, we will give a transparent and intuitive interpretation of what is this ergodic solution and how it is approached in this long time horizon regime.

B.

Overview of recent generalizations.

The class of Mean Field Games described in the previous subsection includes already a rich variety of models which behaviors are far from being fully understood.

In the last few years, very significant effort have been made to extend the Mean Field Games approach beyond this original class of problems, or to consider them from a rather different point of view. These efforts reflects in part the necessity to approach a more realistic description of economic or social (or other) questions, which imposes to relax somehow the restrictive hypothesis originally assumed in the original formulation described above.

As we have already stressed, the behavior of even the simplest Mean Field Games model are very poorly understood (in the sense a physicist uses for this word), and we shall not address in any details these extensions here. Starting from section III, we shall actually further limit our study to a particular subclass of these Mean Field Games (namely the quadratic MFG), which, as we shall see, provide a nice entry point to the field for physicists. We shall however in this subsection and the following one provide a brief survey of the main existing lines of research. We refer the interested reader to recent mathematical review [START_REF] Gomes | Mean field games models -a brief survey[END_REF] and monographs [START_REF] Carmona | Probabilistic Theory of Mean Field Games with Applications[END_REF][START_REF] Carmona | Probabilistic Theory of Mean Field Games with Applications[END_REF].

Generalized cost functions

A very natural way to generalize the basic MFG of section II A is to enlarge the set of admissible cost functions. A first obvious step, actually already taken in the first papers of Lasry and Lions [START_REF] Lasry | Jeux à champ moyen. II -Horizon fini et contrôle optimal[END_REF][START_REF] Lasry | Jeux à champ moyen. I -Le cas stationnaire[END_REF], is to assume that the functional Ṽ in Eq. (2.2) may also have an explicit dependence in time, reading thus Ṽ [m τ ](X i τ , τ ). There are furthermore some situation where the running cost cannot be written as the sum of a "free Lagrangian term" L(xτ, aτ ) depending only on the state and control variable of a given agent and a "potential" term Ṽ [m τ ](x) which describe its interaction with the distribution of all agents. This is the case, for instance when considering pedestrian flow [START_REF] Lachapelle | On a mean field game approach modeling congestion and aversion in pedestrian crowds[END_REF], where congestion effects need to be taken into account. Note that congestion here does not means that agents avoid crowded places (which in any case would be taken into account by a proper choice of the functional Ṽ ), but that being at a location where the density of agents is high is making more costly the use of a large velocity. In that case, it is therefore necessary to introduce a term coupling directly the control parameter with the agent densities.

A particular model for the congestion phenomena where the running cost contains a term proportional to

(m τ + µ) γ |a i τ | β β-1
(with µ, γ ≥ 0 and β ∈ (1, 2]) have been discussed by Lions in [START_REF] Lions | Cours au collège de france[END_REF] and more recently by Adchou and Porretta [START_REF] Achdou | Mean field games with congestion[END_REF].

In another context, Gomes et al. [START_REF] Diogo | On the existence of classical solutions for stationary extended mean field games[END_REF] have also introduced, and analyzed in the stationary case, a model in which the cost function of a given player is affected not only by the distribution of the other players state variable, but also by the value of their control variables. This kind of models arise naturally in the context of Mean Field Game model of trade crowding, and have been coined Mean Field Games of Control by Cardialaguet and Lehalle [START_REF] Cardaliaguet | Mean field game of controls and an application to trade crowding[END_REF].

Different kinds of players

One simplifying assumption in the basic MFG of section II A is that all the player are essentially identical (exchangeable), and therefore distinguish themselves only through the value of their state variable.

A natural extension within the Mean Field game hypothesis is to consider that there exists different groups of agents, G (k) , k = 1, • • • , K, each being characterized by a specific cost function which may also depend separately on the partial densities

m (k) (x) = 1 N i∈G (k) δ(x -X i )
and not only on the total density m(x) = k m (k) (x). In such cases, each agent has to minimize a cost functional similar to Eq. (2.2), but specific to its group, in which

the "potential" Ṽ [m](x) is replaced by Ṽ (k) [m (1) , • • • , m (K) ](x).
One example of these Mean Field Games involving more than one kind of population is the model of segregation studied by Adchou et al. [1], which is in some sense the analogue of the model introduce by Schelling in 1971 to study segregation effects in the United States cities [START_REF] Schelling | Dynamic models of segregation[END_REF].

If the introduction of different kinds of "small players" may lead to rather new kinds of behavior, it however does not lead to a very significant conceptual change in the general theory of Mean Field Games. This aspect is however rather different when the group(s) of small players interact (strategically) with one or many "big players" [START_REF] Carmona | A probabilistic approach to mean field games with major and minor players[END_REF], a situation which is often encountered when modeling financial markets.

For instance, in [START_REF] Lachapelle | Efficiency of the price formation process in presence of high frequency participants: a mean field game analysis[END_REF], Lachapelle et al. considered the question of the price formation in a financial market where a small number of "institutional investors" cohabit with a large number of high frequency traders. In this context, the high frequency traders sell or buy only small quantities and only the coherent action of a large (O(1)) fraction of them can affect the outcome. On the other hand the institutional investors can buy or sell very large quantities that can have a significant impact on their own.

In such circumstance, one can still consider that the mean field hypothesis is still valid for the population of small players. However, the situation is different for the big players as the fully stochastic nature of their individual evolution need to be taken into account. This eventually leads to a mixed description, with a mean field game coupled to the stochastic differential game for a small number of (big) players.

Mean Field Games on Graphs

Another extension of the MFG concept is to apply it to cases where the state space is discrete rather that continuous. This changes rather significantly the structure of the Mean Field Game since the Fokker Planck equation is replaced by a set of rate equations, and similarly the Hamilton-Jacobi-Bellman equation by discrete Bellman equations. The approach to such models is thus technically rather different but may provide a useful simplified setting when analyzing a particular qualitative effect, such as for instance congestion [START_REF] Guéant | Existence and uniqueness result for mean field games with congestion effect on graphs[END_REF]. Furthermore, some economic or social problem naturally lead to such graph-MFG model. In health science for instance, they provide a natural setting for a game theoretical version of the SIR (Susceptible/Infected/Recovered) model for spread of disease, which have been investigated in details by Laguzet and Turinici [START_REF] Laguzet | Individual vaccination as Nash equilibrium in a SIR model with application to the 2009-2010 influenza a (H1N1) epidemic in france[END_REF].

C. Probabilistic approach and the Master Equation

All the Mean Field Game extensions described in the previous subsection involve non-trivial changes with respect to the original formulation of section II A. However, essentially all of them suppose that the evolution of the density of (small) players is described by a deterministic equation and that its fluctuations can be neglected.

However, a number of circumstances require to take into account the existence of stochastic effects, which may survive even in the limit of infinitely many players.

For instance, this is the case in presence of a common noise [START_REF] Carmona | Mean field games with common noise[END_REF] or when one has to take into account the stochastic dynamics of big players [START_REF] Carmona | A probabilistic approach to mean field games with major and minor players[END_REF]. Moreover, such an extension is also needed when one wants to relate an N -player game with its MFG counterpart in the N → ∞ limit.

In these cases, it is not possible to assume that the density of agents can be replaced by its average value, and its stochastic nature needs to be taken into account. At an heuristic level this can be understood as implying that the system of coupled forward Fokker-Planck equation and backward Hamilton-Jacobi-Bellman equation Eqs. (2.5)-(2.6) have to be replaced by a system of coupled stochastic Fokker-Planck equation and stochastic backward Hamilton-Jacobi-Bellman equation. In some instances, such as the case of a large population in a random environment, one can actually construct, and deal with, the resulting stochastic mean field game [START_REF] Carmona | The master equation for large population equilibriums[END_REF]. However, for most cases, it is necessary to reformulate the Mean Field Games approach to take into account the full complexity of these models.

This means that one does not rely any more on an Hamilton Jacobi Bellman equation, which suppose that the time dependent probability distribution of player is fixed, but rather consider this distribution as a variable to be fixed in the optimization process itself. In [START_REF] Carmona | Probabilistic analysis of mean-field games[END_REF], Carmona and Delarue have initiated a purely probabilistic description of Mean Field Games, which is an alternative approach to the one based on PDEs introduced by Lasry and Lions. This probabilistic approach provides one way to tackle the difficulty implied by the stochasticity of the distribution of players [START_REF] Carmona | A probabilistic weak formulation of mean field games and applications[END_REF][START_REF] Carmona | Probabilistic Theory of Mean Field Games with Applications[END_REF]. Another framework for this reformulation has been introduced by P.-L.

Lions and called the "Master Equation" approach [START_REF] Lions | Cours au collège de france[END_REF]. Its basic ingredient consists in writing the value function u (cf Eq. (2.4)) as a function of time, state variable and the full distribution m. This lead to a single differential equation for u which is of second order in m [START_REF] Bensoussan | The master equation in mean field theory[END_REF][START_REF] Cardaliaguet | The master equation and the convergence problem in mean field games[END_REF].

This Master Equation is a very powerful and sophisticated tool, that for obvious reason we shall not describe in any details here, but it represent presumably the most active direction in the development of MFG in the mathematical community [START_REF] Bensoussan | The master equation in mean field theory[END_REF][START_REF] Bensoussan | On the interpretation of the master equation[END_REF][START_REF] Carmona | The master equation for large population equilibriums[END_REF][START_REF] Carmona | Probabilistic Theory of Mean Field Games with Applications[END_REF][START_REF] Carmona | Mean field games with common noise[END_REF][START_REF] Lions | Cours au collège de france[END_REF].

III. QUADRATIC MEAN FIELD GAMES : SCHRÖDINGER APPROACH

A. Quadratic Mean Field Games

In the rest of this paper, we shall restrict our attention to quadratic Mean Field Games, which are defined by the fact that the "free Lagrangian" part of the running cost has a quadratic dependence in the control:

L(X, a) = 1 2 µa 2 (3.1)
This class of MFG will be shown to admit a mapping to a non linear Schrödinger Equation, which leads to an almost complete description of their behavior.

In addition, we will assume that the potential Ṽ [m](x) can be written as the sum of two terms

Ṽ [m](x) = U 0 (x) + V [m](x) (3.2)
where U 0 (x) is an "external potential" which depends only on the state x of the agent while V [m](x) describes the interactions between agents and is invariant under simultaneous translation of both x and m(•). The simplest form that we will consider for this interaction will be linear and local,

V [m](x) = g m(x) (3.3)
where g > 0 corresponds to attractive interactions. We shall also consider non local interactions V [m](x) = dy κ(xy)m(y) , (3.4) and non linear ones

V [m](x) = f [m(x)] , (3.5)
both admitting the simplest form, Eq. (3.3) as a particular case, with, respectively, κ(xy) = g δ(xy) and f (m) = g m.

In the context of crowd dynamics, U 0 (x) would represent the preference of an agent for a given position x, whereas the term V [m](x) takes into account his preference or aversion for crowded places. In this paper we will limit our study to the attractive case (e.g. g ≥ 0 in Equation (3.3)). Two limiting regimes will be of particular interest: the case of strong interactions dominated by V [m](x) and the case of weak interactions in which U 0 (x) is the larger term.

To summarise, we consider a set of N agents, whose individual states at time t are described by continuous variables X i t ∈ R d , which evolves through a controlled linear Langevin dynamics

dX i t = a i t dt + σdW i t , (3.6)
where σ > 0 is a constant, the components of W i are independent white noises of variance 1 and a i is the control chosen by agent i to minimize the cost functional

c[a i ](x i t , t) = T t µ 2 a i τ 2 -Ṽ [m τ ](X i τ ) dτ noise + c T (X i T ) noise (3.7)
where Ṽ [m](x) is a functional of the density m. In this setting, the optimal control is a * t (x) = -1 µ ∇ ∇ ∇u(x, t) with u(x, t) the value function (2.4) and the MFG system (2.5)-(2.6) writes here

∂ t u t (x) - 1 2µ ∇ ∇ ∇u t (x) 2 + σ 2 2 ∆ u t (x) = Ṽ [m t ](x) (3.8) ∂ t m t (x) - 1 µ ∇ ∇ ∇.(m t (x)∇ ∇ ∇u t (x)) - σ 2 2 ∆ m t (x) = 0 (3.9)
with, respectively, final and initial conditions, u T (x) = c T (x) and m 0 (x) = m 0 (x).

In the following we introduce first a change of variables which shows that this system of equations is equivalent to a Schrödinger Equation in imaginary time, and the related formalism which we will use in the rest of this work. We also briefly review three solvable models: non-interacting agents, interactions in the absence of an external potential, and quadratic potential. These models are interesting in their own rights, but they may also serve as reference models in perturbative approaches studied in the next sections.

B. Schrödinger formalism

As a first step, we use the well known fact that the Hamilton-Jacobi-Bellman equation for the value function u(x, t) in (3.8) can be cast into a standard heat equation using a Cole-Hopf transformation [START_REF] Hopf | The partial differential equation u t + uu x = u xx[END_REF] Φ(x, t) = exp -u t (x)/µσ 2 .

(3.10)

The new variable Φ(x, t) obeys a time-backwards diffusion equation,

-µσ 2 ∂ t Φ(x, t) = µσ 4 2 ∆Φ(x, t) + Ṽ [m t ](x)Φ(x, t) , (3.11)
with the final condition Φ(x, T ) = exp(-c T (x)/µσ 2 ). Note that it follows from equation (3.11) that Φ(•, •) > 0 as soon as Φ(x, T ) > 0 everywhere.

The next step is a change of variables for the density m t (x) [START_REF] Guéant | Mean field games equations with quadratic hamiltonian: a specific approach[END_REF] Γ(x, t) = m t (x) Φ(x, t) . (3.12) This second variable now follows a similar heat equation, but forward in time,

µσ 2 ∂ t Γ(x, t) = µσ 4 2 ∆Γ(x, t) + Ṽ [m t ](x)Γ(x, t) , ( 3.13)
with the initial condition Γ(x, 0) = m 0 (x)/φ(x, 0).

Under these transformations, the MFG system has been recast in a pair of nonlinear heat equations, differing only by the sign on the left hand side and by their asymmetric boundary conditions.

Let us now consider the scalar nonlinear Schrödinger equation describing the quantum evolution of a wave amplitude in a reversed potential -Ṽ [ρ],

i ∂ t ψ(x, t) = - 2 2µ ∆ψ(x, t) -Ṽ [ρ](x)ψ(x, t) , ( 3.14)
with ρ ≡ ψ * ψ and R d ρ(x) = 1. We note that Eq. (3.14) and its complex conjugate are equivalent to Eqs. (3.11)-(3.13) under the formal correspondence µσ 2 → ,

φ(x, t) → ψ(x, it) and Γ(x, t) → [ψ(x, it)].
Furthermore the ergodic system (2.7), reads here

-λ e - 1 2µ ∇u e (x) 2 + σ 2 2 ∆u e (x) = Ṽ [m e ](x) (3.15) 1 µ ∇(m e (x)∇u e (x)) + σ 2 2 ∆m e (x) = 0 . (3.16)
When considering the new variables, Φ e = exp (-u e /µσ 2 )) and Γ e = m e /Φ e , we see that both have to follow the same equation

λ e ψ e = - µσ 4 2 ∆ψ e -Ṽ [m e ](x)ψ e , ( 3.17)
with either ψ e = Φ e (x) or ψ e = Γ e . In this context, the connection with the nonlinear Schrödinger equation (3.14) appears more clearly since if ψ e (x) is a solution of Eq. (3.17), then the two time dependent functions

φ(x, t) = exp + 1 µσ 2 λ e t ψ e (x) (3.18) Γ(x, t) = exp - 1 µσ 2 λ e t ψ e (x) , (3.19)
are solutions of, respectively Eqs. (3.11)-(3.13), and simultaneously both solutions of Eq. (3.17), in the very same way as

ψ(x, t) = exp - i λ e t ψ e (x)
is solution of both Eq. (3.14) and Eq. (3.17).

The non linear Schrödinger equation has been used for decades to describe systems of interacting bosons in the mean field approximation (see e.g. [START_REF] Kaup | Perturbation theory for solitons in optical fibers[END_REF][START_REF] Kivshar | Dynamics of solitons in nearly integrable systems[END_REF][START_REF] Kosevich | Particle and wave properties of solitons[END_REF][START_REF] Pérez-García | Dynamics of Bose-Einstein condensates: Variational solutions of the Gross-Pitaevskii equations[END_REF][START_REF] Pitaevskii | Bose-Einstein Condensation[END_REF]), and in the context of fluid mechanics ( [START_REF] Kharif | Rogue Waves in the Ocean[END_REF]). We now introduce a formalism which is well known in these domains and will prove again very useful in the present context of quadratic mean field games.

C. Ehrenfest's relations and conservation laws

By analogy with the NLS Equation [START_REF] Kosevich | Particle and wave properties of solitons[END_REF], we introduce an operator formalism on some appropriate functional space and derive the evolution equations for some quantities of interest. To do so, we first define a position operator X = ( X1 ,

• • • , Xd ),
where Xν acts as a multiplication by the ν th coordinate x ν . We also define a momentum operator Π ≡ -µσ 2 ∇ ∇ ∇. For an arbitrary operator Ô defined in terms of X and Π, we define its average

Ô (t) ≡ Φ(t)| Ô|Γ(t) = dx Φ(x, t) Ô Γ(x, t) , ( 3.20)
where the couple (Φ(t), Γ(t)) defines the state of the system and evolves according to Eqs (3.11)- (3.13). Note that when Ô depends only on the position: Ô = Ô( X), the latter average reduces to the usual mean value with respect to the density:

Ô (t) = dx m t (x)O(x) . (3.21)
Differentiating Equation (3.20) with respect to t, one gets, as for the Schrödinger equation [30], the time evolution of the mean value of an observable in terms of a commutator:

d dt Ô = ∂ Ô ∂t - 1 µσ 2 [ Ô, Ĥ] , (3.22)
where we have introduced the Hamiltonian

Ĥ = - Π2 2µ -Ṽ [m t ]( X) . ( 3.23)
In particular, one gets

d dt X = Π µ , ( 3.24)
d dt Π = F[m t ] , (3.25)
where

F[m t ] ≡ -∇ Ṽ [m t ]( X)
is named by analogy the "force" operator. In the same way, introducing the variance Σ ν of the ν th coordinate, (ν = 1,

• • • , d) , Σ ν ≡ Xν 2 -Xν 2 (3.26)
and the averaged "position-momentum" correlator for the ν th coordinate

Λ ν ≡ Xν Πν + Πν Xν -2 Xν Πν (3.27)
one has

d dt Σ ν = 1 2µ Λ ν Σ ν , ( 3.28)
d dt Λ ν = 2 Xν Fν [m t ] -Xν Fν [m t ] + 2 µ Π2 ν -Π2 ν 2 . (3.29)
Furthermore, when local interactions are assumed,

Ṽ [m](x) = U 0 (x) + f [m(x)] , (3.30)
the mean force depends only on the external potential

F = F0 (3.31) X F = XF 0 -dx xm t (x)f [m t (x)] , (3.32)
where F0 ≡ -∇ x U 0 ( X) , and Eqs. (3.25)-(3.29) can be simplified accordingly.

Finally, with such interactions, we can introduce an action functional

S[Φ, Γ] ≡ T 0 dt R d dx - µσ 2 2 Φ(∂ t Γ) -(∂ t Φ)Γ - µσ 4 2 ∇Φ.∇Γ + Φ U 0 (x) Γ + F [Φ Γ] , (3.33)
where F (m) = m f (m )dm , which extremals are solutions of the system (3.11, 3.13).

Indeed, the variational equation δS/δΓ = 0 (respectively δS/δΦ = 0 ) is equivalent to Eq. (3.11) (respectively Eq. (3.13)) with Ṽ [m] in the form Eq. (3.30). This property will provide us with a variational approximation scheme for the solutions of the MFG system. Note however that the boundary conditions associated with equations (3.11) and (3.13) have to be carefully taken into account. For a pair (Φ, Γ) solving the MFG equations, the action Eq. (3.33) can be rewritten as

S[Φ, Γ] = -dt dx µσ 2 2 Φ(∂ t Γ) -(∂ t Φ)Γ + dt E tot (t)
where we have introduced the "total energy",

E tot (t) ≡ E kin + E pot + E int (3.34)
with "kinetic", "potential", and "interaction" energies defined respectively as

E kin = 1 2µ Π2 , ( 3.35
)

E pot = U 0 (x) , ( 3.36
)

E int = dxF [m(x, t)]. (3.37)
The integrand in the action functional (3.33) does not depend explicitly on time, so that by Noether theorem, there is a conserved quantity along the trajectories. This conserved Noether charge is (due to the sign conventions here) minus the previously defined total energy, so that

dE tot (t) dt = 0 (3.38)
and for any pair (Φ(x, t), Γ(x, t)) solving the MFG equations, one has

S[Φ, Γ] = - µσ 2 2 dt dx [Φ(∂ t Γ) -(∂ t Φ)Γ] + E tot T . ( 3.39)

D. Exactly solvable cases

The list of completely solvable MFG models is up to now rather short, and mainly restricted to stationary settings [START_REF] Gomes | One-dimensional stationary mean-field games with local coupling[END_REF], this situation being most probably due to the difficulties encountered when working within the original representation (2.5)-(2.6).

Hereafter, we add a few examples to that list, by considering first situations in which either the external potential U 0 (x) or the interaction term V [m] is absent or fully negligible. These cases will be also used in the following sections as starting points to develop a perturbative approach when both terms are present in the potential but one is significantly larger than the other and dominates the optimization process.

We conclude the present contribution to that list with the case of an harmonic external potential and local interactions, for which exact solutions can be found for rather specific boundary conditions but no constraint on the relative strength between the two terms of the potential.

Non interacting case

We consider first the uncoupled case V [m] = 0 so that the potential reduces to U 0 (x). The absence of interactions means that the MFG is not really a game anymore since the agents behave independently of the strategies of the others. However, this degenerate case appears naturally in perturbative approaches to weakly interacting regimes which will be considered in Section VI.

Here the potential reduces to the density independent part, Ṽ [m](x) = U 0 (x) so that the Eqs. (3.11) and (3.13) become actually linear,

µσ 2 ∂ t Φ = + Ĥ0 Φ, (Φ(x, T ) = Φ T (x)) (3.40) µσ 2 ∂ t Γ = -Ĥ0 Γ, (Γ(x, 0) = m 0 (x) Φ(x, 0) ) (3.41)
with Ĥ0 the linear operator

H 0 = - Π2 2µ -U 0 (x) . (3.42)
In this case, the solutions of the system Eqs. (3.40)-(3.41) can be expressed in terms of the eigenvalues λ 0 ≤ λ 1 ≤ . . . and the associated eigenvectors ψ 0 (x), ψ 1 (x), . . . of Ĥ0 . Explicitly one has

       Φ(x, t) = ϕ 0 e - λ 0 (T -t) µσ 2 ψ 0 (x) + ϕ 1 e - λ 1 (T -t) µσ 2 ψ 1 (x) + . . . Γ(x, t) = γ 0 e - λ 0 t µσ 2 ψ 0 (x) + γ 1 e - λ 1 t µσ 2 ψ 1 (x) + . . . .
Assuming Ĥ0 non degenerate and the eigenfunctions ψ k normalized, the coefficients

{ϕ k } k=1,2,••• are fixed by boundary conditions at t = T , ϕ k = dx ψ k (x)Φ T (x) , ( 3.43)
which in particular specifies the expression of Φ(x, 0),

Φ(x, 0) = ϕ 0 e - λ 0 T µσ 2 ψ 0 (x) + ϕ 1 e - λ 1 T µσ 2 ψ 1 (x) + . . . . (3.44)
This then fixes the initial value Γ(x, 0), and thus the coefficients

{γ k } k=0,1,••• of Γ(x, t)
as

γ k = dx ψ k (x) m 0 (x) Φ(x, 0) . (3.45)
This spectral analysis allows us to see how, in the non-interacting case, the asymptotic solution converges to the ergodic solution away from the time boundaries when the horizon T becomes very large. Indeed, introducing the characteristic convergence time

τ erg = µσ 2 /(λ 1 -λ 0 ) one has Γ(x, t) γ 0 e - λ 0 t µσ 2 ψ 0 (x) for all t τ erg (3.46) Φ(x, t) ϕ 0 e - λ 0 (T -t) µσ 2 ψ 0 (x) for all t T -τ erg . (3.47)
Hence, when both conditions are fulfilled, the density m(x, t) becomes asymptotically time independent as,

m(x, t) γ 0 ϕ 0 e -λ 0 T σ 2 Ψ 2 0 (x) for all τ erg t T -τ erg .
Normalization imposes that γ 0 ϕ 0 e -λ 0 T σ 2 = 1, so that in the limit of large optimization time, the density profile converges exponentially fast (with the characteristic time τ erg ) to a time independent profile: lim

T →∞ m(x, t) -m e (x) ≤ Ce - t τ erg (3.48) lim T →∞ m(x, T -t) -m e (x) ≤ Ce - t τ erg (3.49)
for all time t, with C a constant. Furthermore, the solution of the ergodic problem Eq. (2.7) for a non interacting mean field game is given by λ e = λ 0 , u e (x) = -µσ 2 log ψ 0 (x) + c and m e (x) = ψ 2 0 (x). Thus any choice of an Hamiltonian Ĥ0 with explicitly known eigenstates would lead to an exactly solvable non-interacting mean field game problem. The list of analytically diagonalisable Ĥ0 indeed contains quite a few systems, among which the case of quadratic potentials. Furthermore, for one-dimensional systems, and more generally for classically integrable Hamiltonian of arbitrary dimensions, very good approximations can be obtained based on the EBK approximation scheme [START_REF] Keller | Corrected Bohr-Sommerfeld quantum conditions for nonseparable systems[END_REF]. We also note here that the ergodic problem requires only the knowledge of the eigenstate ψ 0 associated with the smallest eigenvalue λ 0 , and the rate of convergence to it depends only on the first two eigenvalues.

Local attractive interactions in the absence of external potential.

We now turn to the opposite cases when the external potential U 0 (x) is negligible with respect to interactions. More specifically we consider one dimensional models in which the interaction term Eq. (3.5) is local, with the particular form

V [m](x, t) = f [m(x, t)] = g m(x, t) α
with α > 0 and g > 0. It includes the simple linear form of the interaction potential Eq. (3.3) for α = 1. In such cases, the stationary (ergodic) problem Eq. (3.17) reduces to a generalized Gross-Pitaevskii equation

- µσ 4 2 ∂ 2 xx ψ e -g(ψ e ) 2α+1 = λ e ψ e .
(3.50)

The lowest energy state can be computed using a known procedure [START_REF] Pitaevskii | Bose-Einstein Condensation[END_REF], that we recall for convenience in Appendix B. It is associated with an energy

λ e = - 1 4 Γ( 2 α ) Γ( 1 α ) 2 2α 2-α g α + 1 2 2-α 2α µσ 4 α 2-α (3.51) ( Γ(•)
is the Euler's Gamma function), and has the following expression

ψ e (x) = ψ M cosh( x -x 0 η α ) -1 α , ( 3.52)
where the maximum value Ψ M reads

ψ M = λ e (α + 1) g 1/2α . ( 3.53)
The stationary solution is a localized density around some arbitrary point x 0 (a soliton in the language of the NLS equation), and its typical spatial extension

η α = 2 √ α Γ( 1 α ) 2 Γ( 2 α ) α 2-α α + 1 2α µσ 4 g 1 2-α (3.54)
depends only on the ratio (µσ 4 /g) and results from the competition between the noise which tends to broaden the distribution and the attractive interactions. In the particular case α = 1, the interaction potential becomes linear (Eq. (3.3)) and the above expressions reduce to:

λ e | α=1 = - g 2 8µσ 4 , ( 3.55)
Ψ M | α=1 = g 4µσ 4 ,
(3.56)

η 1 = 2µσ 4 g . (3.57)
Two remarks are in order here. First, the expressions above are clearly not well-defined for α = 2. As we shall discuss in section V, this is related to the fact that the soliton is unstable for α > 2. Moreover, we stress that the generalized Gross-Pitaevskii equation Eq. (3.50) is invariant under translation and therefore the soliton Eq. (3.52) can be centered around any point x 0 of the real axis. In presence of a weak but non zero external potential U 0 (section IV) and local interactions as above, it will follows that a very good approximation for the ergodic state will be a soliton centered at the maximum x max of U 0 . We also use this property hereafter to derive exact results in the case of an external quadratic potential.

Quadratic external potential

In the rest of this section on exact results, we shall use the formal connection

(x) = -1 2 kx 2 , k > 0.
We thus consider a total potential

Ṽ [m](x) = - k 2 x 2 + f (m(x)) .
Following [START_REF] Pitaevskii | Bose-Einstein Condensation[END_REF], we use for Φ and Γ the ansatz

Φ(x, t) = exp - γ(t) -xP (t) µσ 2 ψ e (x -X(t)) (3.58) Γ(x, t) = exp + γ(t) -xP (t) µσ 2 ψ e (x -X(t)) , (3.59)
where ψ e (x) is the solution of the ergodic equation

- µσ 4 2 ∂ 2 xx ψ e (x) -Ṽ [(ψ e ) 2 ](x)ψ e (x) = λ e ψ e (x) ,
which, for small k, is well approximated by the expression on the r.h.s. of Eq. (3.52)

(with x 0 = 0). Note that for this ansatz the resulting density is m(

x, t) = Φ(x, t)Γ(x, t) = ψ 2 e (x -X(t))
, and is thus independent of P (t) and γ(t). Inserting these expressions into the system (3.11-3.13), we get the necessary and sufficient conditions for Eqs. (3.58)-(3.59) to be an exact solution of the time dependent problem:

Ṗ (t) = kX(t) (3.60) Ẋ(t) = P (t) µ (3.61) γ(t) = k 2 X(t) 2 + P (t) 2 2µ -λ e . ( 3.62)
The two first equations describe the motion of the centre of mass X(t) of the density distribution. The third one can be integrated in

γ(t) = X(t)P (t) 2 -λ e t + γ 0 .
This solution describes the evolution of a density distribution with finite spatial extension, that we may call "soliton" because it moves without deformation as a classical particle of mass µ in an inverted quadratic potential U 0 (x) = -1 2 kx 2 . It corresponds, however, to rather specific boundary conditions since the initial density should be of the form m(x, t = 0) = ψ 2 e (x-x 0 ), and the function Φ T specifying the terminal boundary condition should be of the form Φ T (x) = K exp{x p T /µσ 2 }ψ e (x -x T ) where p T and x T are related through the mixed condition x T cosh(ωT ) -p T √ kµ sinh(ωT ) = x 0 .

The two constants γ 0 and K being unessential, this family of solutions is fully described by only two parameters, says x 0 and x T .

Assuming that initial and final conditions have been chosen as above, positions of the center of mass X(t) at initial and final times are then fixed to X(0) = x 0 and X(T ) = x T , and for all intermediate times we get

X(t) = x 0 sinh(ω(T -t)) sinh(ωT ) + x T sinh(ωt) sinh(ωT )
with ω ≡ k/µ. In the long horizon limit T → ∞, apart from initial and final time intervals of order τ erg = 1/ω, the center of mass remains localized in a close vicinity of the unstable fixed point of the external potential U 0 . This is a general feature that we shall discuss in more details in the following section.

IV. STRONGLY ATTRACTIVE SHORT RANGED INTERACTIONS I

This section is devoted to simple one dimensional models where the agents have a strong incentive to coordinate themselves. This first example of asymptotic regime allows us to make a clear exposition of the main concepts that can be effectively used, leading to a rather complete understanding of the behavior of mean field game equations. In the next two sections, we shall use essentially the same tools, addressing somewhat more intricate settings in the same regime in section V, and considering other asymptotic regimes in section VI.

We consider here one dimensional models with interaction potentials which are local and linear as in Equation (3.3). The total potential is therefore of the form

Ṽ [m](x) = U 0 (x) + g m(x) , (g > 0)
, with a weak external potential U 0 (x), in a sense explicited below. We also assume that the initial distribution of agents m 0 (x) is localized and well described by its mean position and variance. We postpone to the next section the discussion on different interaction potentials or initial conditions.

A characteristic feature which can be easily found in the regime of strongly attractive short-ranged interactions is that the agents have a strong incentive to form compact groups evolving coherently which, by analogy with the NLS nomenclature, we shall call "solitons" [START_REF] Martin | Bright solitary-matter-wave collisions in a harmonic trap: Regimes of solitonlike behavior[END_REF]. The initial and final boundary conditions will eventually be an obstruction for the existence of these solitons for a short period of time τ * close to t = 0 and t = T , with τ * → 0 in the limit when the interaction strength g → ∞.

However, for a sufficiently large time horizon T , we expect that the dynamics of such

A. Dynamics of the solitons

In the limit of large interaction strength g → ∞, and excluding a neighborhood of time boundaries, we can assume a strongly localized density of agent m(x, t) with a short characteristic length η. Indeed, for a strength g large enough, the variations of U 0 on the scale η can be considered as weak, and in particular the variations of the external potential around any point, δU 0 = ∇U 0 .δx + γ,γ (∂ 2 γ,γ U 0 )δx γ δx γ + • • • , are dominated by the first term ∇U 0 .δx for a displacement of order |δx| ∼ η. In that case, denoting X t = X (t) the average position of the soliton, and P t = Π (t) its average momentum, the Ehrenfest relations Eqs. (3.24)

d dt P t = F 0 (X t ) -∇U 0 (X t ) . (4.2)
We again recognize the classical dynamics of a point particle of mass µ evolving in the potential U 0 (x) as in the particular example of the quadratic external potential studied at the end of previous section.

However, unlike classical mechanics, where a given trajectory is fully specified by its initial position and momentum, a mean field game problem is defined through mixed initial and terminal conditions. In the present setting, initial formation and final destruction of a soliton should occur fast enough that neither the position of the density center of mass nor the mean momentum are expected to evolve in any significant way in the meanwhile. We are thus led to the following identifications

X t=0 = dx x m 0 (x) (4.3)
P t=T = Π (T ) = µσ 2 dx (∇Φ(x, T )) Γ(x, T ) (4.4)
Eq. (4.3) fixes the initial position of the trajectory; Eq. (4.4) can be written as

P T = -dx m T (x)∇u T (x) (4.5)
which, using the final boundary condition in Eq. (2.5), gives

P T = -∇c T (x) -∇c T (X T ) (4.6)
where the last approximation holds if m(x, T ) is localized on the scale of variations of c T (x), which has to be checked afterwards for consistency.

The dynamics of the soliton is thus the classical dynamics of a point particle of mass µ evolving in the potential U 0 (x), with an initial condition Eq. The mixed initial-terminal character of the boundary conditions have also an implication in the context of the ergodic problem Eq. (2.7) studied in [START_REF] Cardaliaguet | Long time average of mean field games with a nonlocal coupling[END_REF], and to its relationship with unstable fixed points of the dynamics. This is presumably a very general feature of the MFG behavior in the limit of large optimization times T → ∞, and we consider this question in some details here since soliton dynamics is the simplest setting in which it appears.

The dynamics described by Eqs. (4.1) and (4.2) is illustrated in Figure 1 for a one dimensional MFG system for an external potential with a single maximum.

Let τ erg the inverse of the Lyapunov exponent associated with the unstable fixed point {X = X max , P = 0} of the dynamics (X max is defined as the position of the maximum of the potential U 0 and is equal to 0 in this particular example). For large enough values of the time horizon, T τ erg , the MFG system has to spend most of the time in a neighborhood of the fixed point, and the dynamics between initial and final conditions is dominated by the associated stable and unstable manifolds, respectively,

W s = {(X(0), P (0)) ∈ R such that (X(t), P (t)) → (X max , 0) as t → +∞} (4.7) W u = {(X(0), P (0)) ∈ R such that (X(t), P (t)) → (X max , 0) as t → -∞} . (4.8)
In fact, trajectories are essentially identical for all T τ erg : They start from the intersection of the stable manifold W s with the line of initial condition X = X 0 , closely follow W s until they reach a very small neighborhood of the fixed point (X, P ) = (X max , 0) and then switch to the unstable manifold W u that they follow until they reach the intersection of W u with the line P = -∇C T (X) specifying the terminal conditions. Of course for large but finite T the actual trajectories are For large values of T , the soliton has to closely follow the stable manifold (blue curve) up to a small neighborhood of the unstable fixed point at (0, 0) and switch to the unstable one (green line) up to its final position. Here

U 0 (x) = -1 2 x 2 -1 4 x 4 , µ = 1, X 0 = -3/2, c T (x) = x(x -7)/2
slightly off W s and W u and the larger T the closer to these manifolds they are, and thus the longer they stay in the immediate neighborhood of (X max , P = 0). Moreover, the dynamics before and after the dominant fixed point essentially decouple: the dynamics on the stable manifold toward the unstable fixed point is barely affected by a change in the terminal condition at t = T .

For U 0 with a single maxima, the soliton at rest centered on this maxima can be identified with the ergodic state defined in [START_REF] Cardaliaguet | Long time average of mean field games with a nonlocal coupling[END_REF]. When the external potential has multiple local maxima, only the absolute one is associated with the genuine ergodic state, but solitons localized on the other local maxima may play the role of "effective ergodic states" in some configurations.

B. Initial and final stages: formation and destruction of the soliton

We turn now to the slightly more delicate question of describing the formation of the soliton near t = 0 and its destruction near t = T . We limit ourselves here to the simpler case where the initial density m 0 (x) and the final cost function Φ T (x) can reasonably be described through a Gaussian ansatz, and postpone to section V the case of more general configurations.

Variational method for the Gaussian ansatz

One very effective approximation scheme for the Non-Linear Schrödinger equation is the variational method [START_REF] Pérez-García | Dynamics of Bose-Einstein condensates: Variational solutions of the Gross-Pitaevskii equations[END_REF]. A valid approach to it is the use of the action functional Eq. (3.33), from which the system (3.11)-(3.13) can be derived. Variational approximations then amount to minimizing the action only within a small subclass of functions. Assuming that the agents' density is going to contract rapidly around its mean value and then move as a whole toward the optimal position, we consider the following ansatz which is a generalization of the exact solution (3.58)-(3.59) found for quadratic potential U 0 :

Φ(x, t) = exp -γ t + P t • x µσ 2 1 (2πΣ 2 t ) 1 4 exp - (x -X t ) 2 (2Σ t ) 2 (1 - Λ t µσ 2 ) . (4.9) Γ(x, t) = exp +γ t -P t • x µσ 2 1 (2πΣ 2 t ) 1 4 exp - (x -X t ) 2 (2Σ t ) 2 (1 + Λ t µσ 2 ) , ( 4.10)
Within this ansatz, the resulting density m(x, t) = Γ(x, t)Φ(x, t) reads

m(x, t) = 1 2πΣ 2 t exp - (x -X t ) 2 2(Σ t )
U 0 (x) = dx U 0 (x) m(x, t) . ( 4.18)
From now on, we choose γ t = (Λ t /4) so that the last two terms in the last line of Eq. (4.13) cancel. Minimizing the reduced action functional with respect to a variation of the parameters gives first the evolution equations for X t and P t :

Ẋt = P t µ (4.19) Ṗt = -dx ∇U 0 (x) m(x, t) = -∇U 0 t , ( 4.20)
which have the same form as Eqs. (3.24)-(3.25). The r.h.s. of Eq. (4.20) depends both on X and Σ, and generally couples the motion of the center of mass to the shape of the distribution. However, as soon as m(x, t) is sufficiently narrow with respect to the scale given by the inverse curvature of the potential U 0 , the approximation ∇U 0 t ∇U 0 (X t ) holds, and the dynamics for the center of mass (X t , P t ) decouples from the dynamics of (Σ t , Λ t ). In the rest of this section, we consider a situation where that the distribution m(x, t) is sufficiently narrow at all time so that such a decoupling occurs and focus on the dynamics of (Σ t , Λ t ). In particular, the energy of the center of mass, P 2 t /2µ + U 0 (x) , is separately conserved and can be dropped from the expression for the total energy. The evolution equations for the reduced system (Σ t , Λ t ), are thus: The stationary standard deviation Σ * is thus, up to a numerical constant of order one, equal to the width of the exact soliton solution in the absence of an external potential (3.57). The total energy at the fixed point is

Σt = Λ t 2µΣ t (4.21) Λt = Λ 2 t -µ 2 σ 4 2µΣ 2 t + g 2 √ πΣ t . ( 4
Ẽ * tot = 1 8 √ π g Σ * (4.25) (note that Ẽ * tot = 1 2 Ẽ * int = -Ẽ * kin ).
2. Time evolution of the reduced system (Σ t , Λ t ) in the long horizon limit

Here again, (Σ * , Λ * ) is an unstable fixed point for the dynamics, so that, in the long horizon limit, all trajectories will follow the associated stable and unstable manifolds, which fixes the value of the total energy to Ẽ * tot . Using Eq. (4.21) in the expression for the total energy gives then an autonomous equation for the evolution of Σ t on the stable or unstable manifolds

µ 2 Σ2 t - µσ 4 8Σ 2 t + g 4 √ πΣ t = Ẽ * tot , ( 4.26)
which can be set in the form Σt

Σ * = ∓ 1 τ * 1 - Σ * Σ t , ( 4.27)
where the minus (respectively plus) sign in front of the r.h.s. describes the stable (respectively unstable) manifold. The factor τ * is the characteristic time

τ * = µΣ 2 * 2 Ẽ * tot = √ 4π Σ * v g , ( 4.28)
with

v g ≡ g µσ 2 (4.29)
the characteristic velocity associated with the interactions.

Let's us consider first the formation of the soliton. The initial distribution m 0 (x) fixes the initial standard deviation Σ 0 = dx (xx ) 2 m 0 (x). In terms of the reduced variable q t = Σ t /Σ * , Eq. (4.27) can then be integrated as

F (q 0 ) -F (q t ) = t τ * , ( 4.30)
with F (q) = +q + log |1 -q|.

Eq. (4.30) is an exact implicit solution for the motion along the stable manifold of (Σ * , Λ * ) towards the fixed point. It is interesting however to consider various limiting regimes which are derived straightforwardly from the limiting behavior of the function F (q): F (q) -q 2 for q 1 , (

F (q) + log |1 -q| for q 1 , (4.32) F (q) +q for q 1 . (4.33)

We have thus the three following possible behavior depending on the width of the initial distribution:

• If the initial distribution of agent is much narrower than the length scale Σ * characterizing the soliton, the variance Σ 2 t increases linearly at a rate Σ 2 * /τ * until time τ * after which it converges exponentially to Σ * (with the characteristic time τ * ).

• If the initial distribution of agent is already close to Σ * , it converges exponentially to Σ * with the characteristic time τ * .

• If the initial distribution of agent is much wider than the length scale Σ * characterizing the soliton, the standard deviation Σ t decreases linearly at rate Σ * /τ * up to a time t c = (Σ 0 /Σ * )τ * from which is converges exponentially to Σ * (with the characteristic time τ * ).

Considering now the destruction of the soliton near T , and again foreseeing that the density will remain localized on a scale Σ T that can be assumed small, we can write out a simplified form of the terminal condition for Σ 2 . Indeed, the ansatz Eq. (4.9) for Φ implies that u(x, T ) = -µσ 2 log Φ(x, T ), and thus reads

u(x, T ) = γ T + 1 4 µσ 2 log(2πΣ 2 T ) -P T x + µσ 2 -Λ T 4Σ 2 T (x -X T ) 2 + • • •
Assuming the final density localized, we can make a Taylor expansion for the terminal

condition u(x, T ) = c T (x) near X T u(x, T ) c T (X T ) + dc T dx X T • (x -X T ) + d 2 c T dx 2 X T • (x -X T ) 2 2 + • • • .
Identifying term by term the coefficients of both expansions, we recover from the first order term the terminal condition Eq. (4.6) for the center of mass motion, while the second order term gives

d 2 c T dx 2 (X T ) = µσ 2 -Λ T 2Σ 2 T = µ σ 2 2Σ 2 T - ΣT Σ T = µ τ * Σ * Σ T 2 Σ * Σ T -1 , ( 4.34)
where we have used both the evolution equation (4.27) along the unstable manifold and the relation τ * = 2Σ 2 * /σ 2 . In the reduced notations, q T = Σ T /Σ * , the equation for the terminal width reads

q T -2 q 2 T = τ * µ d 2 c T dx 2 (X T ) . (4.35)
In the strong interaction limit, g → ∞, the characteristic time τ * goes to zero and the right hand side of this equation (which is the only term depending on the terminal cost c T (x)) becomes negligible. Thus, in accordance with the approximation scheme used here, one has to take q T 2, so that the final distribution m(x, T ) has an extension of order 2Σ * and thus small, as anticipated. Accordingly, the variance at large times is given by

F (q t ) = F (q T ) - T -t τ * , ( 4.36)
where q T is the smallest solution of Eq. (4.35), and q t , τ * and F (q) are defined as above.

C. Discussion

In this section, we have considered the strong positive interaction regime under three simplifying assumptions : i) the state space is one dimensional; ii) the interaction between the agents is local and linear (V [m](x) = gm(x)); and iii) the initial distribution of agents m 0 (x) is reasonably well described by a Gaussian.

Under these hypothesis, the image that emerges for the evolution of the agent in the state space is extremely simple. The dynamics is divided in three stages : the initial formation of the soliton near t = 0, its propagation, and its final destruction ). In the long horizon limit T → ∞, this motion is furthermore dominated by the maxima of U 0 (x) which correspond to an unstable fixed point of the dynamics where the "soliton" spends most of its time [START_REF] Cardaliaguet | Long time average of mean field games with a nonlocal coupling[END_REF], while the initial (resp. final) motion takes place along the related stable (resp. unstable) manifold. There may also exists intermediate regimes for the long horizon limit where other unstable fixed points, when they exist, may also show up in the dynamics and play the role of "effective ergodic state".

This propagation phase is flanked by significantly shorter initial and final phases where the soliton is respectively formed and destroyed. When the initial distribution of agent m 0 (x) has an extension Σ 0 of the order of the stationary value Σ * or smaller, a soliton forms within a typical time of order τ * (Eq. (4.28)). For initial extensions Σ 0 much larger than Σ * , the time of formation of the soliton is larger, (Σ 0 /Σ * )τ * . In the final phase, the typical extension of the distribution is of order 2Σ 0 for strong interactions and the soliton always disappears in a time of order τ * .

V. STRONGLY ATTRACTIVE SHORT RANGED INTERACTIONS II

In this section, we continue with the study of the strong positive coordination regime, but we relax some of the simplifying assumptions made in section IV concerning the dimensionality of the space, the form of the interactions, and the shape of the initial distribution of agents.

We will show that the dominant phase of the dynamics, namely the propagation of the soliton, is essentially unaffected (or only trivially affected) by these modifications.

Most of our discussions here will concern the formation and destruction phases of the soliton (and in practice we will essentially focus on the former). This section will be divided in three parts. In the first one, we shall extend the variational approach of section IV to higher dimensionality problems and to different form of the interaction between agents. In a second subsection, we shall discuss the "collapse" of the distribution of agents which does occur for nonlinear local interaction or in higher dimension even with linear interactions. Finally, in the last subsection, we shall discuss the formation of the soliton when the initial distribution m 0 (x) has some structure and cannot be just described by its mean and variance.

A. Gaussian ansatz in higher dimensions and nonlinear interactions

In this subsection, we generalize the variational approach of section IV B 1 to the case where the "state space" of the agents is of dimension higher than one and for non-linear interaction between the agents of the form

V [m](x) = gm α (x), so that Ṽ [m](x) = U 0 (x) + g [m(x)] α , (5.1)
with, as in the previous section, g > 0, α > 0 and U 0 (x) assumed non zero but weak.

We first generalize the variational ansatz Eq. (4.10)-(4.9) to

Φ(x, t) = exp -γ t + P t • x µσ 2 d ν=1 1 (2π(Σ ν t ) 2 ) 1/4 exp - (x ν -X ν t ) 2 (2Σ ν t ) 2 (1 - Λ ν t µσ 2 ) (5.2) Γ(x, t) = exp +γ t -P t • x µσ 2 d ν=1 1 (2π(Σ ν t ) 2 ) 1/4 exp - (x ν -X ν t ) 2 (2Σ ν t ) 2 (1 + Λ ν t µσ 2 ) (5.3)
Inserting these expressions into the action Eq. (3.33) (see appendix C), we get an action functional in the variables (X ν t , P ν t ) {ν=1,••• ,d} and (Σ ν t , Λ ν t ) {ν=1,••• ,d} . We consider first the equations of motion for the X ν t and P ν t :

Ẋν t = P ν t µ (5.4) Ṗ ν t = -∂ ν U 0 (x) t -∂ ν U 0 (X t ) (5.5)
which decouples from the dynamics of the Σ ν 's and Λ ν 's when the approximation in (5.5) is assumed. This approach is valid whenever the density of agents is sufficiently narrow with respect to the inverse curvature of the potential, condition which in the strong positive coordination regime will be fulfilled at almost all time (except possibly within a very short time near t = 0 or near t = T ). As in the one dimensional case, the mean position X and mean momentum P hence follow the motion of a classical particle of mass µ in the external potential U 0 (X), with initial and terminal conditions which are the direct generalization of Eqs. (4.3)-(4.6):

X ν t=0 = d d xx ν m 0 (x) (5.6) P ν t=T = -∂ ν c T (x) T -∂ ν c T (X T ) . (5.7)
This motion can be more complex than in the one-d case since the conservation of the total energy is not sufficient to make the dynamics integrable anymore. However, in the long horizon limit, it is still dominated by the maxima of the potential U 0 (x max )

and takes place very close to the stable and unstable manifolds of the unstable fixed point (X = x max , P = 0).

We assume from now on that the motion of the center of mass decouples from the evolution of the shape of the density profile and consider the dynamics of (Σ ν t , Λ ν t ) {ν=1,••• ,d} alone. We consider the case of a local interaction of the form V [m] = gm α and from the variation of the reduced action, we get the following evolution equations:

Σν = Λ ν 2µΣ ν , (5.8) Λν t = (Λ ν t ) 2 -µ 2 σ 4 2µ(Σ ν t ) 2 + 2gα α + 1 d ν =1 1 √ α + 1(2π) α/2 1 Σ ν t α .
(5.9)

These equations admit a single stationary state (Σ * , Λ * ), with Λ ν * = 0 (5.10)

Σ ν * = Σ * =   4α α + 1 1 (α + 1)(2π) α d/2 g µσ 4   -1/(2-αd)
.

(5.11)

As in section III D 2 for the one dimensional case d = 1, there is no solution for α = 2/d. We understand here this critical value of α as the transition between the situation 0 < α < 2/d where (Σ * , Λ * ) is a saddle point for the total energy of the reduced system Eq. (C.10), and the situation α > 2/d where it is a minima, leading to a change of stability. From a physical point of view, this means that for α > 2/d, attractive interactions dominate at short distance while diffusion, which tends to disperse the density m(x, t) dominates at large distance, which makes the "soliton" unstable. Note that the stability (resp. instability) of the soliton is associated with instability (resp. stability) of trajectories.

B. Collapse for dα > 2

To get a better picture of the main differences between the regime α < 2/d where the soliton is stable and the regime α > 2/d where it is unstable, we restrict ourselves to the 1-dimensional case (thus α c = 2) and introduce the canonical variables

q t = Σ t Σ * , ( 5.12
)

p t = - Σ * 2 Λ t Σ t , ( 5.13)
The Lagrangian Eq. (C.8) reads L(t) = p t qt -h(p t , q t ) , (5.14)

h(p, q) = - p 2 2µΣ 2 * + µσ 4 4Σ 2 * 1 2q 2 - 1 αq α , (5.15)
and the equation of motions takes the canonical form in term of the Hamiltonian

h(p, q) q = + ∂h(p, q) ∂p = - p µΣ 2 * (5.16) ṗ = - ∂h(p, q) ∂q = µσ 4 4Σ 2 * 1 q 3 - 1 q (α+1) .
(5.17)

With these variables, conservation of the total energy Ẽtot = -h(p, q) is manifest, and the fact that the Liouville measure dp dq is conserved (which would be also true for d > 1 or in the full problem when variables (p, q) are coupled with the global motion (P, X)) makes it possible to classify a priory the fixed points by their stability.

Specifically here, the dynamical system Eqs. (5.16)-( 5.17) has one fixed point at (q * = 1, p * = 0), where the second derivatives of h(q, p) are given by

∂ 2 h ∂ 2 p (q * ,p * ) = -1 µΣ 2 * , ∂ 2 h ∂p∂q (q * ,p * ) = 0 , ∂ 2 h ∂ 2 q (q * ,p * ) = µσ 4 4Σ 2 *
(2 -α) . (5.18) We see therefore that, as expected, the stability of the fixed point is entirely determined by the sign of (α -2).

• For 0 < α < 2, (q * , p * ) is a saddle point for h(q, p) and thus an unstable fixed point. The dynamics is in that case qualitatively similar to the one of section IV B 2: Near t = 0 (formation of the soliton), the system starts from the initial q 0 fixed by the initial density m 0 (x) and follow the stable manifold of (q * , p * ), which it approaches exponentially closely on a very short time scale.

The destruction of the soliton follows a similar scenario, but on the unstable manifold. The typical phase portrait in this case is shown in Figure 2a.

• If α > 2, (q * , p * ) is a minima of h(q, p) and thus a stable (elliptic) fixed point for the classical dynamics governed by Eqs. (5.16)-(5.17). For a given set of initial and final conditions, one cannot exclude the possibility that a periodic orbits in the neighborhood of (q * , p * ) turns out to be solution of the equations of motion, which would correspond to a kind of breathing of the soliton. Our guess, however, is that these breathing mode, when they exist, are not the only solution of the equations of motion, and can be eliminated because they do not minimize the cost Eq. (3.7) (in the sense that they correspond to local minima, but not absolute minima of this cost). In the limit of long time horizon, the system will prefer to flow toward other fixed points: either a low density (non-Gaussian) noise-dominated phase (described in our ansatz by the limit q → ∞), or a large density phase dominated by the interactions (here obtained in the limit q → 0). This case is illustrated in Figures 2b-2c in the particular case α = 3.

Hence, for α > 2, we have two possible options: either a large spreading of the distribution, or a collapse. In the first case, namely a large excursion toward large q's (and thus large Σ's) the initial spreading of the density should be large enough so that the noise becomes the dominating force. Within the approximation scheme we use here, what the system does is then to spread out relatively slowly under the influence of the noise, and (possibly) re-compactify toward the end (i.e. for t near T ) if the terminal constraint makes this mandatory. In practice however, a system in this configuration is effectively not any more in the strong interaction regime. There is no short time scale associated with the interaction between the agents, and the influence of U 0 (x) may become as significant as the one of the noise. Furthermore since the distribution of agents does not remain localized, the dynamics of (q t , p t ), does not decouple from the center of mass (X t , P t ), and even the validity of the Gaussian ansatz Eqs. (4.10)-(4.9) becomes questionable. The analysis of this regime should actually follow the line of section VI.

Second, if the initial and final conditions select a regime where the density of agents remains sufficiently large, then the system will rather choose a large excursion toward Σ t → 0, and thus a collapse of the density of agents. In that case, we need to consider explicitly a "finite-range" interaction. Indeed, the rational behind the utilization of a "zero-range" interaction potential Eq. (3.5) is that the actual range of the interactions is the smallest length scales in the problem, which cannot hold any more here since Σ t would eventually become smaller than whatever this range is.

Let us illustrate this in the case α = 3 depicted in Figures 2b-2c. The interaction 3 can be seen as the ξ → 0 limit of

V [m](x) = gm(x)
V [m](x) = g dy 2 dy 3 dy 4 K(x, y 2 , y 3 , y 4 )m(y 2 )m(y 3 )m(y 4 ) , with K(y 1 , y 2 , y 3 , y 4 ) ≡ 1 2( √ 2πξ) 3 exp   - 1 16ξ 2 i =j (y i -y j ) 2   -→ ξ→0 δ(y 1 -y 2 )δ(y 2 -y 3 )δ(y 3 -y 4 ) .
The analysis of this "finite range" interaction can be done along the same line as before, up to the replacement of the interaction energy term by

Ẽint = g 8   1 2π(ξ 2 + Σ 2 t )   3 -→ ξ→0 g 8(2π) 3/2 Σ 3 t
(which is indeed the second term of Eq. (C.10) for d = 1 and α = 3).

With the (q t , p t ) variable, the Hamiltonian Eq. (5.15) becomes

h(p, q) = - p 2 2µΣ 2 * + µσ 4 4Σ 2 *   1 2q 2 - 1 3 1 + ξ 2 Σ 2 * ξ 2 + Σ 2 * ξ 2 + Σ 2 * q 2 3/2  
(which as ξ → 0 indeed correspond to the Hamiltonian Eq. (5.15) with α = 3). Here Σ * is the value of Σ at a stationary point and is a solution of

ξΣ 4 * (ξ 2 + Σ 2 * ) 5/2 = 2(2π) 3/2 3 ξµσ 4 g . ( 5.19)
Note that the left hand side of Eq. (5.19) depends only on the ratio Σ * /ξ and has a single maximum at Σ * /ξ = 2. For ξ (and the right hand side of (5.19)) small enough, Eq. (5.19) has thus exactly two positive solutions, says, Σ

(1) * and Σ

(2) * , which are respectively smaller and larger than 2ξ, each one associated to a stationary point for the dynamics; the second derivative of H(p, q) with respect to q now reads

∂ 2 h ∂ 2 q (q * ,p * ) = µσ 4 4Σ 2 * 4ξ 2 -Σ 2 * ξ 2 + Σ 2 * (5.20)
while the other two second derivatives remain as in (5.18). Thus the smallest value Σ

(1) * is associated with an hyperbolic fixed point and the larger one Σ

(2) * with an elliptic fixed point. The corresponding phase portrait is shown on Figure 3, where both fixed points appear. The stable "soliton" is associated with this new fixed point, and its size Σ

(1) * is governed by the range of the interaction ξ and not any more by the balance between the strength of the interaction and the one of the noise.

C. Non-Gaussian initial densities of agents

In this last part of the section, we shall consider the situation where we relax the assumption that the initial density of agents m 0 (x) can be correctly described by a Gaussian. For sake of clarity we limit this discussion to the one-dimensional case and local linear interactions Ṽ [m](x) = +gm(x). We first consider the situation where the initial density m 0 (x) can be described as the juxtaposition of two well separated Gaussian-like bumps, and then discuss some aspects of the general case. Furthermore, we restrict the discussion to the initial times (formation of the soliton), assuming that final boundary conditions are compatible with a Gaussian distribution.

To clarify the question of structured, but non Gaussian, initial distributions of agents, it is useful to consider the example of an initial condition which can be split into two well separated parts m 0 (x) = m a 0 (x) + m b 0 (x), both separately well approximated by a Gaussian and characterized by their relative masses, mean positions and standard deviations, hereafter denoted by ρ k = m k 0 (x)dx, X k 0 and Σ k 0 , with k = a, b, respectively. We consider a variational ansatz which is a straightforward generalization of Eqs. (4.10)-(4.9), namely

Φ(x, t) = Φ a (x, t) + Φ b (x, t) , Γ(x, t) = Γ a (x, t) + Γ b (x, t) ,
where, for k = a, b,

Φ k (x, t) = ρ k exp -γ t + P k t • x µσ 2 1 2π(Σ k t ) 2 1/4 exp - (x -X k t ) 2 (2Σ k t ) 2 (1 - Λ k t µσ 2 ) Γ k (x, t) = ρ k exp +γ t -P k t • x µσ 2 1 2π(Σ k t ) 2 1/4 exp - (x -X k t ) 2 (2Σ k t ) 2 (1 + Λ k t µσ 2 )
and follow the same approach as in section IV B 1. We make the two following assumptions: i) both parts m a t (x) and m b t (x) remain well separated for the time necessary to get an equilibrium shape, |X a t -X b t | Σ a t + Σ b t , and ii) the extensions Σ a t , Σ b t of both parts are small on the scale at which U 0 (x) varies significantly. Under these assumptions, the time evolution greatly simplifies as each sub-part behaves independently form the other and follow Eqs. (4.1)-(4.2) for its center of mass, and Eqs. (4.21)-(4.22) for the distribution parameters, with an effective coupling constant g k = ρ k g in the interaction potential. Thus, using the results of the previous section, we get that each sub-part k forms a soliton with a rescaled extension Σ k * = 1 ρ k Σ * (with Σ * as in (4.23)), which implies that the smaller part gets the larger extension. This first evolution takes place on time scales τ k * = ρ k τ * (see Eq. (4.28) and the discussion at the end of section IV B 2). Let us now analyze the separate motion of the centers of mass for each part, (X k t , P k t ). In order to keep with the simplest picture, we completely neglect the effect of U 0 (x) (relaxing this assumption does not introduce major conceptual changes) and suppose that in the long time limit, the density of agents should form a single soliton of mass one and at rest. In such a case, the solitons evolve as independent classical particles with unknown constant velocities v k = P k 0 /µ, which have to be determined.

forms a local soliton of extension inversely proportional to its mass which moves until it merges with a neighboring soliton. Furthermore, if the size ρ k of these sub-parts is uniformly bounded from below, then the formation of local solitons is characterized by the velocity scale v g and occurs on a short time scale in the limit of strong interactions. In this setting, lighter solitons take more time to form, and move faster than the heavier ones. When more than two solitons are present, various scenarios are possible which differ by the order in which they merge together, implying different choices of initial conditions.

More general initial densities of agents with inhomogeneities but no clearly separated sub-parts would have to be studied in a case by case basis. However, the fact that the extension of local solitons is inversely proportional to their mass lead us to expect the formation of solitons for strong enough interaction potential, even if the determination of their distribution would remain a-priori a difficult problem.

VI. PERTURBATIVE APPROACH TO THE WEAKLY INTERACTING

REGIME

In contrast with the two previous sections, we now turn to the case when the interactions between agents are small with respect to the external potential so that they can be described as a perturbation of a non-interacting model.

The general strategy here is relatively clear. If the "interaction potential"

V [m](x)
is small, one should first solve the (non-interacting) Schrödinger equation (as in section III D 1); we then plug in the interactions, assumed to be small, and insert the potential term V [m](x) as a perturbation, using the standard tools of quantum mechanics. We shall see however that the forward/backward structure of the Mean Field Game equations introduces some subtleties in this relatively straightforward scenario.

Here, we limit ourselves to the description of interactions up to first order corrections. Furthermore, we shall consider here the long horizon limit T → ∞, and concentrate mainly on the convergence to the ergodic state. We thus discard the effects of the final boundary conditions, which show up only in the late stage of the process.

As an application, we will consider a one-dimensional model with a quadratic (inverted) external potential

U 0 (x) = - k 2 x 2 (6.1)
and a weak short-ranged interaction potential

V [m](x) = g m(x) (6.2)
where g is a small positive coupling constant.

A. Non-interacting model

We start with a brief discussion of the non-interacting limit Ṽ [m](x) ≡ U 0 (x), mainly to fix some notations and to recall some properties we shall make use of in this section.

From the results of section III D 1, we know that the time evolution of both functions Φ and Γ can be derived from the eigenfunctions ψ n (x) and eigenvalues λ n of the Hamiltonian

H 0 = - 1 2µ Π2 -U 0 (x) . (6.
3)

The time evolution of the two functions Φ(x, t) and Γ(x, t) can be expressed in terms of these eigenfunctions through the construction of a propagator

G H 0 (x, x , t) ≡ n≥0 e -λn t/µσ 2 ψ n (x)ψ n (x ) , ( 6.4)
where the subscript H 0 stands for the "free" Hamiltonian (6.3). We have

Φ(x, t) = dx Φ(x , t ) G H 0 (x , x, t -t) t ≤ t (6.5) Γ(x, t) = dx G H 0 (x, x , t -t ) Γ(x, t ) t ≥ t (6.6)
We now consider the influence of both the initial density of agents m 0 (x) = Φ(x, 0)Γ(x, 0) and the terminal condition Φ(x, T ) = K exp (-c T (x)/µσ 2 ). In the long horizon limit T τ erg , where we define the ergodic time as Γ H 0 e (x, t) ≡ C -1 e -λ 0 t/µσ 2 ψ 0 (x) (6.9)

τ erg ≡ µσ 2 λ 1 -λ 0 , ( 6
(with C some arbitrary constant that we fix to one here), so that the resulting density profile is time independent, m(x, t) = m H 0 e (x) = ψ 2 0 (x). The backward time evolution of Φ(x, t) coming from the ergodic state Φ H 0 e at some fixed final time is trivial, and in particular Φ(x, 0) = ψ 0 (x). The initial condition for Γ thus reads

Γ(x, 0) = m 0 (x) ψ 0 (x) (6.10)
and for all further times t

Γ(x, t) = dx G H 0 (x, x , t) m 0 (x ) ψ 0 (x ) . ( 6.11)
Since m(x, t) = Φ(x, t)Γ(x, t) (cf (3.12)), the time evolution for the density can be written as m(x, t) = dx F H 0 (x, x , t)m 0 (x ) (6.12)

where we have introduced the density time-propagator F H 0 (x, x , t) ≡ ψ 0 (x)G H 0 (x, x , t) e +λ 0 t/µσ 2 ψ 0 (x ) . (6.13)

As stressed before, these expressions for the propagation of the density of agents are valid in the long optimization time limit T τ erg , and their simplicity can be eventually traced back to the fact that Φ(x, t) remains in its ergodic state Φ e as long as (T -t) τ erg , and in particular near t = 0.

We end this subsection with a few comments. First, we stress that for times large enough, the propagator (6.4) factorizes G H 0 (x, x , t) e -λ 0 t/µσ 2 ψ 0 (x)ψ 0 (x ) for all t τ erg , (6.14) and one recovers as expected the ergodic state Γ(x, t) = Γ e (x, t) from (6.11), and m(x, t) = m H e (x) from (6.12). This implies in particular that for any normalized initial density m 0 (x)

dx F H 0 (x, x , t)m 0 (x ) m e (x)
for all t τ erg . (6.15)

Finally, one can check easily that m e is a fixed point of the propagation equation (6.12), m H 0 e (x) = dx F H 0 (x, x , t)m H 0 e (x ) , (6.16) as again expected. We shall make use of these properties below.

B. First order perturbations : the ergodic state

We want now to compute the first order corrections to the previous non-interacting model when a weak interaction term V (x, t) = g m(x, t) is added to the potential.

Let us denote by ψ e the solution of the nonlinear ergodic problem, λ e ψ e (x) = -µσ 4 2 ∆ψ e (x) -U 0 (x)ψ e (x) -g (ψ e (x)) 3 . (6.17)

For g small enough, both ψ e and λ e can be computed using perturbation theory around the lowest energy state of H 0 (Eq. ( 6.3)). To first order in g, one gets easily

ψ e (x) = ψ 0 (x) + g n>0 V 0,n (λ n -λ 0 ) ψ n (x) + o(g) (6.18
)

λ e = λ 0 -g V 0,0 + o(g) (6.19)
where ψ n is the unperturbed eigenfunction of H 0 and for all n, m ≥ 0,

V n,m ≡ dx ψ 2 0 (x)ψ n (x)ψ m (x) . (6.20)
The density in the ergodic state then reads, up to first order,

m e (x) = ψ 2 0 (x) + 2g ψ 0 (x) n>0 V 0n (λ n -λ 0 )
ψ n (x) + o(g) .

C. First order perturbations : dynamics

We now construct the dynamic evolution toward the ergodic state, given an initial density profile m 0 .

The basic tool we shall use in this subsection is essentially the time dependent perturbation theory of quantum mechanics. However the forward/backward structure of the mean field game equations introduces some extra complication since, as we shall see, it requires to have perturbative results which remain valid for very long times (of the order of the optimization time T ).

For this reason, it turns out to be necessary to develop the perturbation theory not around the unperturbed Hamiltonian H 0 (Eq. (6.3)) but around the Hamiltonian associated with the true ergodic state

H e (x) = - Π 2 2µ -U 0 (x) -g (ψ e (x)) 2 , ( 6.21)
where for now ψ e (x) denotes the exact solution of the ergodic problem Eq. (6.17).

We then write the full time dependent Hamiltonian as H e plus a perturbation:

H(x, t) = H e (x) -g (m(x, t) -m e (x)) (6.22)
where m e (x) is the (assumed known) ergodic density of state and m(x, t) is the yet unknown time dependent density of agents.

Let T e be a time at which the system is in the ergodic state (possibly up to an exponentially small error in T that we fully neglect here). Φ(x, t) is solution of Note that in this last expression, the reasons for the choice of a perturbation theory around H e rather than around H 0 can be made clear: first, the exponential factor may differ greatly from the same expression with λ 0 instead of λ e , since t can get large independently of g; second, with the present choice, the time integral in the right hand side is well defined, even in the limit T e → ∞ (assuming the convergence of the expansion).

µσ 2 ∂ t Φ(x, t) = H(x, t)Φ(x,
In particular, the value at t = 0 reads Φ(x, 0) = 1 + g µσ where we have used also the fact that F e (x, y, t) m e (x) to get the last line, valid up to first order in g. For short times, t ≤ µσ 2 /|V 1,1 -V 0,0 | (assuming that V 1,1 = V 0,0 ), the dynamics towards the ergodic state does not differ too much from the unperturbed one and the densities and propagators in the right hand side of Equation (6.30) can be replaced by their first order approximations in the first line, and their expressions at g = 0 in the next three lines. We thus get an explicit form the first order solution to the Mean Field Game equations:

m(x, t) = m H 0 e (x) + dx F H 0 (x, x , t) (m 0 (x ) -m H 0 e ) + dx (F e (x, x , t) -F H 0 (x, x , t)) m 0 (x ) + g µσ 2 t 0 ds dy dx [F H 0 (x, y, t -s) -F H 0 (x, x , t)] × m H 0 (y, s) -m H 0 e (y) F H 0 (y, x , s)m 0 (x ) + g µσ 2 Te t ds dy dx m H 0 (y, s) -m H 0 e (y) × [F H 0 (y, x, s -t) -F H 0 (y, x , s )] F H 0 (x, x , t)m 0 (x ) (6.32)
where the term on the second line accounts for the first order correction to the propagator.

For large times, (t τ erg ), this expression converges exponentially to m e (x) for the same reasons as in Eq. (6.31).

For short times, (t

τ erg ), it leads to m(x, t) = m 0 (x) + ∂ t m(x, 0) t + O ((t/τ erg ) 2 ),
where, within first order approximation,

∂ t m(x, 0) ∂ t m (0) 0 (x) + g µσ 2 ∂ t m (1) 0 (x) (6.33) with ∂ t m (0) 0 (x) = dx ∂ t F H 0 (x, x , 0)m 0 (x ) (6.34)
the "free" contribution. We can furthermore write the first order correction as the sum

∂ t m (1) 0 (x) = ∂ t m (e) 0 (x) + ∂ t m (a) 0 (x) + ∂ t m (r) 0 (x) (6.35)
with

∂ t m (e) 0 (x) = g µσ 2 -1 dx [∂ t F e (x, x , 0) -∂ t F H 0 (x, x , 0)]m 0 (x ) (6.36) ∂ t m (a) 0 (x) = -m 0 (x)
∂ t F H 0 (x, y, 0) = σ 2 2 (ψ 0 (x)δ"(x -y) -ψ" 0 (x)δ(x -y)) 1 ψ 0 (y) ,
(and the equivalent expression for ∂ t F e (x, y, 0)).

The term ∂ t m (e) 0 can be understood as originating from the modification of the ergodic state (cf section VI B), ∂ t m (a) 0 is related to the influence of interactions on anticipations and ∂ t m (r) 0 is due to the more (in time) classical, retarded effect of interactions.

We now apply these results to the example of the harmonic oscillator potential Eq. (6.1), for which we compute the first order, O(g/µσ 2 ), corrections to ∂ t m(x, 0).

D. Weakly interacting agents in an harmonic potential

We consider now in more details the harmonic case

U 0 (x) = - k 2 x 2 .
Our goal here is to obtain explicit expressions for the various quantities involved in Eq. (6.32), and more specifically the density propagator F H 0 and the time-dependent density of agents m H 0 (x, t).

When the potential is harmonic the eigenfunctions of the unperturbed Hamiltonian (6.3) can be written as

ψ n (x) = 1 √ 2 n n! 1 π 1/4 √ 0 exp - x 2 2 2 0 H n x 0 , n ≥ 0 (6.40) where o = σ µ k 1/4
and H n (u) is the n th Hermite polynomial; the associated eigenvalues are λ n = µσ 2 ω(n + 1/2) with ω = k/µ. In particular the ground state of H 0 reads

ψ 0 (x) = 1 π 1/4 1/2 o exp - 1 2 x 2 2 o . ( 6.41)
By Mehler formula (see reference [START_REF] Watson | Notes on generating functions of polynomials: (2) Hermite polynomials[END_REF] or Appendix D for a derivation) the propagator Eq. (6.4) reads explicitly

G H 0 (x, x , t) = 1 2π 2 0 sinh(ωt) exp - (x 2 + x 2 ) cosh(ωt) -2xx 2 2 0 sinh(ωt) (6.42) e -ωt/2 √ πl 0 exp - (x 2 + x 2 )) 2 2 0 when (ωt 1)
This in turns implies that

F H 0 (x, x , t) = G Σ F (t) (x -x e -ωt ) , (6.43) with Σ F (t) = 0 (1 -e -2ωt )/2 , where G Σ is a centered Gaussian of width Σ G Σ (x) = 1 √ 2πΣ exp - x 2 2Σ 2 . ( 6.44)
Using Eq. (6.43), the implementation of Eq. (6.32), beyond the zero'th order term, now reduces to quadratures.

In the particular case of a Gaussian initial condition

m 0 (x) ≡ G Σ 0 (x -x 0 ) (6.45)
the integration in Eq. (6.12) can be performed and the time dependent density profile is Gaussian at all times, with

m H 0 (x, t) = G Σm(t) (x -x 0 e -ωt ) (6.46) with Σ m (t) = 0 (1 -(1 -2 Σ 2 0 / 2 0 ) e -2ωt ) /2 . (6.47)
For time small enough, the density-propagator F H 0 (x, x , t) is peaked around x = xe ωt , and the integral in Eq. (6.12) is dominated by a neighborhood of size e ωt Σ F (t) around this value. If the initial profile m 0 (x ) is slowly varying on this length-scale, the corresponding term can be factorized out of the integral in Eq. (6.12) (which is akin to performing a stationary phase approximation), leading to a simpler expression for the density at short times m(x, t) e ωt m 0 xe ωt .

(Note that this does not require m 0 to be a Gaussian.) If the typical scale of variations Σ 0 of the initial distribution of agents is significantly larger than the length 0 characterizing the ground state ψ 0 (x), this approximation will be valid up to time t ∼ ω -1 log(Σ 0 / 0 ).

We turn now to the short time evolution of the density for a Gaussian initial density as in (6.45). In expression (6.33) all integrations over space are Gaussian and can be made explicitly, which leaves only the integration over time to be performed numerically. The result up to first order in the interaction strength is illustrated in Fig. 4 for a particular set of parameters (see caption for more details).

A few remarks are in order. In Fig. 4a, the short time free evolution of the density is what is expected, i.e. a motion toward the maxima x = 0 of the "utility" function U 0 (x); the corrections due to the interactions are indeed quite small in the case presented here. This smallness is due in part to the fact that the three contributions in (6.35) "push" in different directions and compensate one another for a large part (see Fig. 4c). The terms ∂ t m state and ∂ t m (a) 0 (x) associated with anticipations tend to accelerate the motion, when the retarded contribution ∂ t m (r) 0 (x) tends to slow it down (as the interactions make the initial mean location slightly more favorable than in the free case). Beyond a the way.

In this paper, we have mostly limited ourselves to an introduction of these methods in the context of Mean Field Games, and illustrated them on a few simple examples, but it is already clear that this is very far from exhausting the possible applications of the connection with the nonlinear Schrödinger equation.

To start with, preliminary results show that other tools developed in the context of the nonlinear Schrödinger equations, from simple ones such as the Thomas Fermi approximation to more sophisticated one related to inverse scattering methods, can be used to analyze further the behavior of "paradigmatic" population models in the same spirit as what have been done here. More elaborated models, including for instance the presence of two or more groups of agents with different behaviors, or other specific modifications can also certainly be analyzed following the same approach.

Obviously, the class of models considered here forms only a small subset of all possible Mean Field Games, but is already large enough to contain non trivial examples for which no explicit exact solutions are to be found. It includes both potential models (for which an action functional can be defined) and non-potential ones (for which this is not possible), as well as monotone (repulsive interactions in the model treated here ) and non-monotone ones (here the attractive case) [59].

Therefore, it seems to us relatively clear that, within a relatively short amount of time, the connection between quadratic Mean Field Games and the nonlinear Schrödinger equation will provide a good understanding of a large class of Mean Field Game problems. We believe that quadratic MFG are at some level representative of a much more general class of Mean Field Games and expect that the connection outlined here will contribute significantly to a better understanding of general MFG models, as well as extending the field of their possible applications. In this process, we are confident that the physics community, with its specific knowledge and point of view, can, and should, play an significant role. We hope that this paper will be instrumental in this respect. which can be readily integrated as

Integrating once gives

ψ e (x) = ψ M cosh( x -x 0 η α ) -1 α . (B.7)
Finally, the value of λ is fixed through the normalization condition:

dx (ψ e (x)) 2 = 1 .

Setting

I α = +∞ 0 (cosh(x)) -2 α dx = 4 1-α α Γ( 1 α ) 2 Γ( 2 α ) , (B.8)
we get

λ e = - g α + 1 2 2-α α √ 2µσ 4 I α 2α 2-α = - 1 4 Γ( 2 α ) Γ( 1 α ) 2

  with the NLS equation(3.14) to derive particular exact solutions of the MFG system (3.11)-(3.13) in dimension one, for a local interaction potential of the form Eq. (3.5) and a quadratic external potential U 0

  solitons dominates for a large time interval of order [τ * , T -τ * ]. This naturally raises a few questions that split in two sets: On the one hand, we have to understand what are the shape and characteristic scales of these solitons, and how and how fast they form near t = 0 and disappear near t = T . On the other hand, and maybe more importantly since it dominates most of the time interval, we need to understand what govern their dynamics. We first address this simpler question on the dynamics of the solitons, and will consider in a second stage their formation and destruction near the time boundaries.

  (4.3) for the position at t = 0 and a mixed terminal condition Eq. (4.6) involving position and momentum at t = T .It should be stressed however that, compared to the classical situation for which the initial position and momentum are specified, such boundary conditions change drastically the qualitative behavior of the system under study. To start with, while the specification of both initial position and momentum entirely determines a trajectory, a finite number of trajectories may fulfill the mixed conditions Eqs. (4.3)-(4.6). One may therefore have to evaluate the cost functional Eq. (3.7) on each of them to select the correct solution of the MFG problem. Furthermore, such a mode of selection indicates that a MFG system may switch abruptly from one type of trajectory to another under a small variation of some parameter and possibly of the optimization time, which would correspond to a genuine phase transition in the MFG behavior.

Figure 1 .

 1 Figure 1. Phase portrait in the plane (X, P ) for the dynamics of the center of mass of a one dimensional MFG soliton in an external potential with a single maximum. The vertical dashed red line is the loci of the initial states compatible with initial position X 0 = -3/2; The slanted dashed line is the loci of final states compatible with the mixed final condition P T + X T = 7/2. The evolution time T determines the actual trajectory.

  near t = T . Actually, for most of the time interval ]0, T [ the density of agents m(x, t) is well approximated by a soliton of extension Σ * ∼ η (1) (Eq. (3.57)), which is the shortest length scale of the problem. This soliton evolves as a classical particle in a potential U 0 (x) (i.e. following the Hamilton equations Eqs. (4.1)-(4.2)), with the initial and terminal conditions Eqs. (4.3)-(4.4

Figure 2 .

 2 Figure 2. Phase portrait in the canonical variables (q, p). (a): 0 < α < 2. The fixed point is unstable and long time trajectories stay close to the stable and unstable manifold (α = 1, µσ 2 = 1 ). (b) and (c): α > 2. The fixed point is locally stable (b) but non closed trajectories can be seen at larger scale (c) (α = 3, µσ 2 = 1).

Figure 3 .

 3 Figure 3. Phase portrait in the canonical variables (q, p) for α = 3 and a non local interaction kernel of range ξ = .3Σ * . Besides the elliptic fixed point at q = 1 (blue dot) already present in the limit ξ → 0, there appears a new (hyperbolic) fixed point at q = O(ξ) (red dot) which governs the dynamics for optimization time long enough (µσ 2 = 1).

Figure 4 .

 4 Figure 4. Short time evolution of a density of agent in an harmonic potential in the limit of weak interactions (the units of time and length are the one set by the harmonics oscillator, i.e. ω = 1 and l 0 = 1). (a) Initial density of agents and its short time evolution. Solid black: Initial density of agent m 0 (x), chosen here as a Gaussian of width Σ = .5 centered at x 0 = 2. Solid gray: Short time evolution m 0 (x) + ∆t∂ t m (0) 0 (x) in the non interacting approximation (∆t = 0.1). Dashed black: Short time evolution m 0 (x) + ∆t(∂ t m (0) 0 (x) + (g/µσ 2 )∂ t m (1) 0 (x)) including the first order corrections associated with interactions ( the "small parameter" (g/µσ 2 ) has been set to one to enhance visibility). (b) Time derivative of the density. Solid gray: non interacting term ∂ t m (0) 0 (x). Solid black first order interaction term ∂ t m (1) 0 (x). (c) The three contributions to the interacting term (Eq. (6.35)). Dashed (blue online) contribution ∂ t m (e) 0 (x)) associated with the perturbation of the ergodic state. Long-dashed (magenta online) contribution ∂ t m (a) 0 (x)) associated with anticipations. Dotted (brown

µσ 4 4 (Defining a characteristic length η α as η α = µσ 4 - 2 η 2 α∂ x ψ e (x) ψ M 2 -

 44222 ∂ x ψ e (x)) 2 + g 2α + 2 (ψ e (x)) (2α+2) + λ 2 (ψ e (x)) 2 = 0 , (B.2)where the integration constant (the right hand side of Equation (B.2)) is set to zero since a solution associated with the minimum value for λ e has to decay to zero at infinity: lim |x|→∞ ψ e (x) = lim |x|→∞ ∂ x ψ e (x) = 0 . The function ψ e (x) need to have (at least) a nonzero maximum ψ M , which value is thus the unique positive solution of g 2α 2 λ e , (B.5) equation (B.2) can be reduced in the form α

  the expression of the characteristic length (B.5) becomesη α = 2 Γ( 1 α ) 2 Γ( 2 α) external potential as in (5.1).

7 ) 10 )

 710 The two pairs of variables, (X, P) and (Σ, Λ) are coupled only through the potential energy (C.6), as a consequence of the curvature of U 0 on the scale Σ t .Assuming that these corrections are negligible, the two pairs decouple and evolve independently. The motion of the center of mass follows a reduced dynamics in the external potential U 0 :Ẋt = P t µ , Ṗt = -∇U 0 (X) ,and the total energy of the center of mass P 2 2µ + U 0 (X) is separately conserved. On the other hand, the dynamics in the center of mass, (Σ t , Λ t ) is governed by the reduced action S(Σ, Λ) = T 0 dt L(t) withL(t) = Lτ (t) + Ẽtot (t) , t ) 2 -µ 2 σ 4 8µ(Σ ν t )The 2d equations of motion (5.8)-(5.9) are obtained by computing the variations of the action along the trajectories and equating them to zero. one stationary point at which the variances and positionmomentum correlators are the same in all directions, namely Λ ν * = Λ * = 0 energy at the stationary point (Σ * , Λ * Λ t and its derivatives from the evolution equations (C.11)-(C.12) lead to a set of second order coupled equations in the variables Σ ν t only. (C.14)-(C.15), and introducing the reduced variables q ν t = Σ ν t /Σ * , we get the simpler expression:

  and (3.25) together with

	Eq. (3.31) reduce to
	d dt	X t =	P t µ	,	(4.1)

  Inserting that variational ansatz into the action Eq. (3.33), we get S =

							T	L(t)dt
	with the Lagrangian L = Lτ + Ẽtot where		0
	Lτ = -dx	µσ 2 2	[Φ(x, t)(∂ t Γ(x, t)) -(∂ t Φ(x, t))Γ(x, t)]	(4.13)
	= Ṗt X t -	Λ t 2Σ t	Σt -γt +	1 4	Λt	(4.14)
	with the [conserved] total energy
		Ẽtot ≡ Ẽkin + Ẽint + U 0 (x) ,	(4.15)
	where
		Ẽkin =	P 2 t 2µ	+	Λ 2 t -µ 2 σ 4 t 8µΣ 2	,	(4.16)
		Ẽint =	g √ 4 πΣ t	.	(4.17)
							2	,
	which is a Gaussian centered in X t with standard deviation Σ t . Furthermore, for Γ
	and Φ given by Eqs. (4.10)-(4.9),
				Π = P t ,	(4.11)
				Λ = Λ t ,	(4.12)

with Λ defined by Eq. (3.27). P t is thus the average momentum at time t, and Λ t the average position-momentum correlator of the system.

  .7) the system gets close to the ergodic state at all intermediate times t such that t τ erg and (T -t) τ erg . The terminal condition becomes thus irrelevant, except possibly in the late stages that we do not consider here. The ergodic state in the absence of

	interactions is
	Φ H 0 e (x, t) ≡ C e +λ 0 t/µσ 2 ψ 0 (x)	(6.8)

  (y) e (λes/µσ 2 ) [m(y, s) -m e (y)]G e (y, x, s -t)where we have denoted by G e (y, x, t) the propagator (6.4) associated with H e .

			t)	(6.23)
	with the terminal condition
		Φ(x, T e ) = e λeTe/µσ 2 ψ e (x) .	(6.24)
	Following standard time-dependent quantum perturbation theory [63], Φ(x, t) reads,
	up to first order in the perturbation and for all t in [0, T e ],
	Φ(x, t) = e λet/µσ 2 ψ e (x)
	+ ds dy ψ e In analogy with Eq. (6.13), we introduce the density time-propagator associated Te g µσ 2 t
	with H e ,
	F e (y, x, t) = ψ e (y) G e (y, x, t)e +λet/µσ 2 1 ψ e (x)	,	(6.25)
	and we write the evolution of Φ as
	Φ(x, t) = e λet/µσ 2 1 +	µσ 2 g	Te
				(6.26)

t ds dy [m(y, s) -m e (y)]F e (y, x, s -t) ψ e (x) .

  respectively, the third one is specific to the forward/backward structure and corresponds to a modification of the initial data Γ(•, 0).

	As a coherence check of the above expression, we can consider the long time
	limit t → ∞, of Equation (6.30) and verify that, at first order in g, it is indeed
	coherent with the expected result that m(x, t) m e (x) whenever t τ erg . Using
	the fact that at the lowest order, |m(x, t) -m e (x)| ≤ Ce -t/τerg , (which we write as
	m(x, t) m e (x)), we get from Eq. (6.14), at first order in g
	2 dx F e (x, x , t)m 0 (x ) Te ds dy [m(y, s) -m e (y)]F e (y, x, s) ψ e (x) 0 Given an initial distribution of agents m 0 , the initial value of Γ(x, 0) can be now (6.27) m(x, t) computed up to first order as Γ(x, 0) = 1 -g µσ 2 Te 0 ds dy [m(y, s) -m e (y)]F e (y, x, s) m 0 (x) ψ e (x) (6.28) Now, since Γ(x, t) is solution of -µσ 2 ∂ t Γ(x, t) = H(x, t)Γ(x, t) (6.29) + t g µσ 2 0 ds dy dx F -t g µσ 2 0 ds dydx F m e (x) (6.31)
	it can be written up to first order in g as
	Γ(x, t) = dx G e (x, x , t)Γ(x , 0)
	+	g µσ 2			0	t	ds dx dy G e (x, y, t -s)[m(y, s) -m e (y)]G e (y, x , s)Γ(x , 0) .
	Collecting the previous results, one can write an expression for the density of
	agents at first order in perturbation theory. We get
	m(x, t) = dx F e (x, x , t) m 0 (x )
	+	g µσ 2	0	t	ds dy dx F e (x, y, t -s) [m(y, s) -m e (y)] F e (y, x , s)m 0 (x )
	+	g µσ 2	t	Te	ds dy [m(y, s) -m e (y)] F e (y, x, s -t)	dx F e (x, x , t)m 0 (x )
	-	g µσ 2	0	Te	ds dydx F e (x, x , t) [m(y, s) -m e (y)] F e (y, x , s )m 0 (x ) (6.30)
	Note that there are thus three terms at the first order of perturbations with a

different origin: though the two last terms terms are rather classical and correspond to the first order perturbation of each of the two factors Γ(x, t) and Φ(x, t), e (x, y, t -s) [m(y, s) -m e (y)] F e (y, x , s)m 0 (x ) e (x, x , t) [m(y, s) -m e (y)] F e (y, x , s )m 0 (x )

dx F e (x, x , t)m 0 (x ) + g µσ 2 t 0 ds dy dx [F e (x, y, t -s) -F e (x, x , t)]

× [m(y, s) -m e (y)] F e (y, x , s)m 0 (x )

  Lagrangian L(t) = L τ + (E kin + E pot + E int ) whereThe computation of the potential energy would requires the external potential U 0 to be given. However, since E pot depends only on m(x, t) = Γ(x, t)Φ(x, t), and not separately on both Γ(x, t) and Φ(x, t) it depends on X t and Σ t only. Furthermore, using a Taylor expansion of U 0 (x) around X t , we get E pot (t) = U 0 (X t ) + 1 2

	and we choose from now on
				γ t = γ 0 +	d ν=1	Λ ν t 4	,	(C.5)
	so that the first two terms in (C.3) cancel.
	The action Eq. (3.33) can be written as
									T
			S[Φ, Γ] =	dtL(t) ,
							d	0
	ν=1 (Σ ν t ) 2 ∂ 2 νν U 0 (X) + O( Σ t Finally, the interaction energy in the variational ansatz reads with a L τ (t) = -µσ 2 2 dxΦ(x, t) (∂ E kin (t) = 1 2µ Π2 t = -µσ 4 dx∇Φ(x, t) • ∇Γ(x, t) , E int (t) = g [Φ(x, t)Γ(x, t)] α+1 dx α + 1 2 d ν=1    1 2π(Σ ν t ) 2 1/4 exp -(x ν -X ν 4 ) . t ) 2 (2Σ ν t ) 2 (1 -E µσ 2 = g α + 1 d α=1   1 2π(Σ ν t ) 2 (α + 1)(x -X ν t ) 2 dx exp 2(Σ ν t ) 2	Λ ν t µσ 2 )	(C.6)    .
									(C.1)
	Γ(x, t) = exp	+γ t -P t • x µσ 2	d ν=1	  	1 2π(Σ ν t ) 2 1/4 exp -	(x ν -X ν t ) 2 (2Σ ν t ) 2 (1 +	Λ ν t µσ 2 )	   ,
									(C.2)
	Computation of the Lagrangian
	We obtain for the first two terms
		L d ν=1	Λ ν t 2Σ ν Σν t t	(C.3)
		E kin (t) =	P 2 t 2µ	+	d ν=1	(Λ ν t ) 2 -µ 2 σ 4 8µ(Σ ν t ) 2 ,	(C.4)

t Γ(x, t)) -(∂ t Φ(x, t))Γ(x, t) , pot (t) = U 0 t = dx Φ(x, t) U 0 (x) Γ(x, t) , E int (t) = dx F [Φ(x, t)Γ(x, t)] = g α + 1 [Φ(x, t)Γ(x, t)] α+1 dx .

We hereafter develop a variational ansatz by minimizing this action on a restricted class of functions (Φ(x, t), Γ(x, t)) defined as in (5.3)-(5.2) that we recall here for convenience:

Φ(x, t) = exp -γ t + P t • x τ (t) = -γt + d ν=1 Λν t 4 + Ṗt • X t -

Acknowledgments

: This research has been conducted as part of the project Labex MME-DII (Funded by the Agence Nationale pour la Recherche, Grant No. ANR11-LBX-0023-01).

Conservation of total energy and the final condition chosen implies that the system of the two initially separated solitons have the same energy as a single soliton at rest, that is

where we have used that the energy of a soliton in the center of mass Ẽ * tot (see Eq. (4.25)) scales as the square of the coupling constant g. Then, in the absence of an external potential U 0 , total momentum is conserved (cf Eq. (3.25)), k=a,b µ ρ k v k = 0 .

(5.22)

The velocities of the solitons before collision are thus

for k = a, b and with v g defined by Eq. (4.29), the velocity scale associated with the interactions.

If the pair of solitons have an extension initially larger than their invariant value Σ k * , they contract with the initial velocity given by Eq. (4.27),

We find that light solitons have a contraction velocity slower than their center of mass, resulting in a positive velocity of the front of matter toward the other soliton before they reach their equilibrium shape. However, equilibration time is smaller for lighter solitons, τ k c = ρ k Σ 0 Σ * τ * , so that the front of matter, X k t + Σ k t moves by a finite fraction of Σ 0 in the time necessary for Σ k t to reach the value Σ k * (the maximal value over ρ k is found to be

3) Σ 0 .27Σ 0 ). Thus the picture we gave here is consistent with the hypothesis of an initial separation, Σ a 0 + Σ b 0 |X a 0 -X b 0 |. For arbitrary initial conditions, the exact scenario may become significantly more complex, and a precise description which would be universally valid is obviously beyond the scope of the present work. We limit ourselves to what can be anticipated on a general basis.

The case of two solitons studied above can easily be generalized to a larger number: if the initial density of agents m 0 (x) can be separated in a few non-overlapping subpart of mass ρ k and size larger than (Σ * /ρ k ), each of these sub-parts contracts and tendency to make the distribution slightly more narrow, the net effect of interactions for this particular example is thus to effectively slightly slow the motion of the group.

We stress however that for the example considered here the width of the initial density has been chosen slightly smaller than the one of the ergodic state, and their mean positions not too far away. Other choices may have led to a stronger influence of the anticipations.

VII. CONCLUSION

In this paper, after a general introduction to Mean Field Games in the form in which they have been introduced orginally by Lasry and Lions, [START_REF] Lasry | Jeux à champ moyen. II -Horizon fini et contrôle optimal[END_REF][START_REF] Lasry | Jeux à champ moyen. I -Le cas stationnaire[END_REF][START_REF] Lasry | Mean field games[END_REF], and a bird's eye survey of the recent mathematical development in that field, we have considered in details a class of mean field game models, referred to here as "quadratic" MFG.

Such models describe the collective behavior of a large number of agents, whose individual dynamics follow a controlled linear Langevin equation, when the control derives from the minimization of a quadratic cost functional.

As we have emphasized, there exist a formal, but deep, relationship between the MFG equations describing these models and the nonlinear Schrödinger equation. Our main purpose was to explore this relationship and its implications on the structure of the solutions of MFG problems. Indeed, the nonlinear Schrödinger equation has a very long history in physics, and many tools and approximation schemes have been developed along the years to analyze its properties in different parameter regimes. Using the connection between MFG and NLS, we have shown that it is possible to adapt some standard tools from Quantum Mechanics (Ehrenfest relations, perturbative expansions), or to rely on concepts and techniques more specific to the nonlinear Schrödinger equation (here the concept of soliton, the existence of an action from which the equations of motion derive, and the related variational approaches) and that it may indeed lead to a sharp control over the behavior of quadratic MFG models in various regimes. In particular we have introduced variational methods that are well adapted in the limit of strong interaction, while on the opposite, weak interaction limit, a perturbation theory can be developed. A few (partly new) exact results have been derived along The Hamilton-Jacobi-Bellman equation is one of the basic tool of optimal control and is discussed in details in many textbooks, such as for instance [START_REF] Bertsekas | Dynamic programming and optimal control[END_REF]. As it might be less familiar to physicists, we provide a brief sketch of its derivation below.

Starting from the definition Eq. (2.4) of the value function u(x, t), the principle of dynamic programming consists in splitting the optimization in two parts, the first one for the infinitesimal time interval [t, t + dt], and the second for all times beyond this. This reads

with dX given by the Langevin equation (2.1). Using Itô lemma we then have

which, inserted in Eq. (A.1) yields the Hamilton-Jacobi-Bellman equation

where H(x, p) ≡ inf α α α (L(x, α α α) + p.α α α)).

The boundary condition u(x, t = T ) = c T (x) is then obtained by noting that at the end of the optimization interval, there is no control variable on which the agent can act to optimize its cost, and therefore the utility function is just given by the final cost c T .

This completely solve the optimization problem for the cost function when one assumes that the value function is twice differential. In more general settings, weak or viscosity solutions for the HJB equation has to be considered instead [START_REF] Lions | Optimal control of stochastic integrals and hamiltonjacobi-bellman equations[END_REF].

Appendix B: Solutions of the generalized Gross-Pitaevskii equation in one dimension

In one dimension and in the absence of external potential, U 0 = 0, the ergodic problem considered in subsection III D 2 reduces to the following generalized Gross-Pitaevskii equation (3.50):

Appendix C: Gaussian variational ansatz

In this appendix, we provide some of the intermediate results that has been used when discussing the variational approach used in sections IV and V.

Invariant manifolds in one dimension

The fixed point (Σ * , Λ * ) of the dynamical system Eqs. (5.8)-(5.9) is unstable for α ∈]0, 2[ and the associated soliton is stable. Along the stable and unstable manifolds associated with the stationary point, the total energy Ẽtot = Ẽ * tot and Equation (C.20) reads

In one space dimension, the expressions for the width (C.14) and the total energy (C.15) at the stationary point simplify:

and the equation (C.21) can be integrated, giving the equation for the stable and unstable manifolds. Introducing the function F (q) as the integral

if q > 1 .

(C.24)

The equation for the stable manifold reads

while the equation for the unstable manifold is

where the characteristic time τ * is

Note that the function F (q) behaves as F (q) ≈ 1 √ α log |1 -q| when q ≈ 1, as F (q) ≈ 2-α 2α q 2 for q 1 and like F (q) ≈ q for q 1 .

Appendix D: Propagator for the harmonic oscillator

In this appendix, we give a brief derivation of the expression for the propagator Eq. (6.42) corresponding to an harmonic potential. This formula is a rather classical result and various derivations can be found in [START_REF] Watson | Notes on generating functions of polynomials: (2) Hermite polynomials[END_REF]; the following one is given here for completeness. We first set the length unit so that 0 ≡ 1.

We consider the n th eigenfunction of the harmonic oscillator ψ n (6.40) and introduce its Wigner transform as

We get

where Reinserting by homogeneity the dependence on the length scale 0 leads to Eq. (6.42).

[1] Y. Achdou, M. Bardi, and Cirant M. Mean field games models of segregation.

Mathematical Models and Methods in Applied Sciences, 27(1):75-113, 2017.