∎ 11institutetext: F. Bagagiolo 22institutetext: Department of Mathematics,
Università di Trento, Italy 22email<EMAIL_ADDRESS>33institutetext:
S. Faggian 44institutetext: Department of Economy, Università Ca’ Foscari
Venezia, Italy 44email<EMAIL_ADDRESS>55institutetext: R. Maggistro
66institutetext: Department of Mathematical Sciences, Politecnico di Torino,
Italy 66email<EMAIL_ADDRESS>77institutetext: R. Pesenti
88institutetext: Department of Management, Università Ca’ Foscari Venezia,
Italy 88email<EMAIL_ADDRESS>
# Optimal control of the mean field game equilibrium for a pedestrian
tourists’ flow model

Fabio Bagagiolo  Silvia Faggian  Rosario Maggistro  Raffaele Pesenti

(Received: date / Accepted: date)

###### Abstract

Art heritage cities are popular tourist destinations but for many of them
overcrowding is becoming an issue. In this paper, we address the problem of
modeling and analytically studying the flow of tourists along the narrow
alleys of the historic center of a heritage city. We initially present a mean
field game model, where both continuous and switching decisional variables are
introduced to respectively describe the position of a tourist and the point of
interest that it may visit. We prove the existence of a mean field game
equilibrium. A mean field game equilibrium is Nash-type equilibrium in the
case of infinitely many players. Then, we study an optimization problem for an
external controller who aims to induce a suitable mean field game equilibrium.

###### Keywords:

Tourist flow optimal control mean field games switching variables dynamics on
networks

###### MSC:

91A13 49L20 90B20 91A80

##  1 Introduction

In the recent years, art heritage cities have experienced a continuous growth
of tourists to the point that overcrowding is becoming an issue and local
authorities start implementing countermeasures. In this paper, we address the
problem of modeling and analytically studying the flow of tourists (more
precisely pedestrian daily excursionists) along the narrow alleys of the
historic center of a heritage city. Starting from the results in Bagagiolo-
Pesenti bagpes , we recast this into a mean field game for controlled dynamics
on a network, representing the possible paths in the city together with the
possible attractive sites.

We suppose that the tourists have only two main attractions to visit. This is
for simplicity only: situations with more attractions can be equally treated.
We represent their possible paths inside the city by a circular network where
three points are identified: $S$ the train station where tourists arrive in
the morning and to which they have to return in the evening; $P_{1}$ the first
attraction; $P_{2}$ the second attraction (see Figure 1.a). We are given an
external arrival flow at the station, represented by a continuous function
$g:[0,T]\to[0,+\infty[$ (roughly speaking, $g(t)$ stays for the density of
arriving tourists per unit time). Here $T>0$ is the final horizon, the time
before which everyone has to be returned to the station, after the tour of the
city.

Figure 1: a) the physical network of the paths inside the city (with the three
identified points: train station, $S$, and two attractions, $P_{1},P2$) and
two tourists visiting the city following opposite directions. b) the state
network and the states $(\theta,w_{1},w_{2})$ of the tourists during their
visits. Each singular circle-branch represents the network of the city paths
as seen by tourist with given values of the switching variables $w_{1},w_{2}$.
The dashed arrows between the four circles-branches represent the switching of
the four labels: at the beginning the label is $(1,1)$; when the attraction
$P_{1}$ is reached (and hence visited), the label switches to $(0,1)$, and
similarly when $P_{2}$ is reached; the last is $(0,0)$ which holds when both
attractions are reached. The point-dashed arrow represents the external
arrival flow in the station.

Each single tourist (agent) controls its own dynamics represented by the
equation

$\theta^{\prime}=u$ (1)

where $\theta\in\mathbb{R}$ is a space-coordinate in the network and $t\mapsto
u(t)\in\mathbb{R}$ is a measurable control. We denote by $\theta_{S}$,
$\theta_{1}$ and $\theta_{2}$, respectively, the position of the station, of
the attraction $P_{1}$ and of the attraction $P_{2}$. To each tourist, we
associate a time-varying label $(w_{1},w_{2})\in\\{0,1\\}\times\\{0,1\\}$. For
$i\in\\{1,2\\}$, $w_{i}(t)=1$ means that, at the time $t$, the attraction
$P_{i}$ is not visited yet, and $w_{i}(t)=0$ that the attraction was already
visited. The state of an agent is then represented by the triple
$(\theta,w_{1},w_{2})$, where $\theta$ is a time-continuous variable and
$w_{1},w_{2}$ are switching variables. In the following, we denote by
$B=[0,2\pi]\times\\{0,1\\}\times\\{0,1\\}$ the state space of variables
$(\theta,w_{1},w_{2})$ and in particular, we call (circle-)branch any subset
$B_{\hat{w}_{1},\hat{w}_{2}}$ of $B$ which includes the states
$(\theta,{\hat{w}_{1},\hat{w}_{2}})$, with $\theta\in[0,2\pi].$ Such branches
also correspond to the edge of the switching networks in Figure 2, which is
another way to represent our network model with its switching dynamics.

Indeed, while the evolution of the value of $\theta$ is governed by (1), the
evolution of the values switching variables can be only from $1$ to $0$ and is
described by the following condition

$w_{i}=\begin{cases}1&\text{for}\ t\in[t,\tau_{i}]\\\ 0&\text{for}\
t\in]\tau_{i},T]\end{cases}$

being $\tau_{i}\in[t,T]$ is the first time instant at which the agent reaches
and visits attraction $P_{i}$, $i\in\\{1,2\\}$. This dynamics is represented
by the arrows and the labels in Figure 1.b (see also Figure 2). In Figure 1.b,
consider a tourist that arrives in the station and visits $P_{1}$ first, at
time $\tau$. It has state $(\theta_{S},1,1)\in B_{1,1}$ when it arrives at the
station. Then, at time $\tau\in[t,T]$, its state $w_{1}$ switches, it was
equal to 1 in $[t,\tau_{1}]$, it assumes value $w_{1}=0$ in a left-open
interval $]\tau_{1},T]$. Immediately after $\tau_{1}$, the tourist’s state
switches to the branch $B_{0,1}$. See also Figure 2.

The cost to be minimized by every agent, when it starts its tour from the
station at time $t$, is given by a combination of terms representing: i) the
hassle of running during the tour; ii) the pain of being entrapped in highly
congested paths; iii) the frustration of not being able to visit some
attractions; iv) the disappointment of not being able to reach the station by
the final horizon $T$. Such a cost can be analytically represented by

$\begin{array}[]{ll}J(t,u)=\int_{t}^{T}\left(\frac{u(s)^{2}}{2}+{\cal
F}^{w_{1}(s),w_{2}(s)}({\cal M}(s))\right)ds+\\\
c_{1}w_{1}(T)+c_{2}w_{2}(T)+c_{S}\xi_{\theta=\theta_{S}}(T)\end{array}$ (2)

Here, $c_{1},c_{2},c_{S}>0$ are fixed, and
$\xi_{\theta=\theta_{S}}(t)\in\\{0,1\\}$ and it is equal to $0$ if and only if
$\theta(t)=\theta_{S}$. In (2), the quadratic term inside the integral stays
for the cost i); the second term inside the integral stays for the congestion
cost ii); the first two addenda outside the integral stay for the cost iii);
the third addendum outside the integral stays for the cost iv). In particular,
the congestion cost ${\cal F}^{w_{1}(s),w_{2}(s)}({\cal M}(s))$ is
instantaneously paid by an agent whose switching label at time $s$ is
$(w_{1}(s),w_{2}(s))$ being the actual distribution of the agents ${\cal
M}(s)$. For any $(w_{1},w_{2})\in\\{0,1\\}\times\\{0,1\\}$, ${\cal
F}^{w_{1},w_{2}}$ is a positive function defined on the set of all admissible
distribution of agents.

In Bagagiolo-Pesenti bagpes  the problem of minimizing the cost (2) by a large
number of tourists is seen as a mean field game and, under some suitable
assumptions, the existence of a mean field game equilibrium is proven. A mean
field game equilibrium is a time-varying distribution of agents on the
network, $t\mapsto{\cal M}^{*}(t)$ for $t\in[0,T]$, such that, when inserted
it in (2), the returned optimal control $u^{*}(\cdot;t):[t,T]\to\mathbb{R}$,
$s\mapsto u(s;t)$, which is implemented by all agents who starts moving from
the station at time $t\in[0,T]$, gives rise to optimal trajectories of the
agents which exactly generate the time-varying distribution ${\cal M}^{*}$. A
mean field game equilibrium can be seen as a fixed point, over a suitable set
of time-varying distributions, of a map of the form

${\cal M}\longrightarrow u_{\cal M}\longrightarrow{\cal M}_{u_{\cal M}}$ (3)

where $u_{\cal M}$ is the optimal control when $\cal M$ is inserted in (2),
and ${\cal M}_{u_{\cal M}}$ is the corresponding evolution of the agents’
distribution when all of them are moving implementing $u_{\cal M}$ as control.
Of course, the problem must be coupled with an initial condition for the
distribution $\cal M$, whereas the boundary condition is represented by the
incoming external arrival flow $g$.

We remark that the concept of mean field equilibrium is a Nash-type
equilibrium concept. Indeed, in the case of a large number of agents, even
infinitely many, as the case of mean field games is, every single agent is
irrelevant, the single agent has measure zero, it is lost in the crowd. Hence,
in the case of equilibrium, for a single agent is not convenient to
unilaterally change behavior, because such a single choice will not change the
mean field $\cal M$, and the agent will not optimize.

The goal of the present paper is twofold. First we amend some stringent
assumptions that were made in bagpes  in order to prove the existence of a
mean field game equilibrium. Then we study a possible optimization problem for
an external controller who aims to induce a suitable mean field game
equilibrium. We suppose that such an external controller (the city
administration, for example) may act on the congestion functions ${\cal
F}^{w_{1},w_{2}}$, choosing them among a suitable set of admissible functions.

Consider as an example the historical center of the city of Venice, Italy.
Tourists typically enter the city at the train station and from there they
have two main alternative routes, a shorter and a longer one, to the main
monuments. In the recent years, the shorter route have been particularly
congested on peak days. For this reason, local authorities have introduced
both some gates to slow down the access to the shorter route and some street
signs to divert at least part of the tourist flow long the longer route.

The mean field games theory goes back to the seminal work by Lasry-Lions
LLions  (see also Huang-Caines-Malhamé HCAINMAL ), where the new standard
terminology of Mean Field Games was introduced. This theory includes methods
and techniques to study differential games with a large population of rational
players and it is based on the assumption that the population influences the
individuals’ strategies through mean field parameters. Several application
domains such as economic, physics, biology and network engineering accommodate
mean field game theoretical models (see Achdou-Camilli-Capuzzo Dolcetta,
AcCamDolc , Lachapelle-Salomon-TuriniciLachSalTur ). In particular, models to
the study of dynamics on networks and/or pedestrian movement can be found for
example in Camilli-Carlini-Marchi CCMar , in Cristiani-Priuli-Tosin CPTos ,
Camilli-De Maio-Tosin CDeMTos  and Bagagiolo-Bauso-Maggistro-Zoppello BBMZop .
Our problem is also a sort of routing problem. Indeed to control the roadway
congestion, different strategies was proposed as variable speed limits (Hegyi-
Schutter-Hellendoorn HegyiSchutter1 ), ramp metering (Gomes-Horowitz GomesHoro
) or signal control (Brian Park-Yun-AhnBrianYun ). However, such mechanisms do
not consider neither the agents’ perspective nor affect the total amount of
vehicles/people. A significant research effort was done to understand the
agents’ answer to external communications from intelligent traveler
information devices (see Srinivasan-Mahmassani Mahamassani , Khattak-
Polydoropoulou-Ben AkivaKhattak ) and, in particular, to study the effect of
such technologies on the agents’ route choice behaviour and on the dynamical
properties of the transportation network (see Como-Savla-Acemoglu-Dahleh-
Frazzoli ComoSavla ). Moreover, it is known that if individual agents make
their own routing decisions to minimize their own experienced delays, overall
network congestion can be considerably higher than if a central planner had
the ability to explicitly direct traffic. From that, the idea to include in
our problem an external controller which induces a suitable mean field game
equilibrium. While, in the specific case of vehicular congestion, tolls
payment is considered to influencing drivers to make routing choices that
result in globally optimal routing, namely to force the Wardrop equilibrium to
align with the system optimum network flow (see Smith Smith , Morrison
Morrison , Dial Dial  and Cole-Dodis-Roughgarden Cole ). The Wardop
equilibrium Wardrop  is a configuration in which the perceived cost associated
to any source-destination path chosen by a nonzero fraction of drivers does
not exceed the perceived cost associated to any other path. This is a
stationary concept, indeed recently was developed its continuous counterpart
(see Carlier-Jimenez-Santambrogio Carlier1 , Carlier-Santambrogio CarlierSant
) which fits the situation of pedestrian congestion (that in this paper we
indeed treat using a mean field model) and also it is useful to look at large
scale traffic problems, when one only want to identify average value of the
traffic congestion in different zones of a large area. Actually, also the
models using that continuous framework are essentially stationary because they
only accounts for sort a cyclical movement, where every path is constantly
occupied by the same density of vehicles, since those who arrive at
destination are immediately replaced by others. This is the essentially
difference with respect to the mean field models. Indeed in the latter, due to
the explicit presence of time, the optimal evolution is given by a system
coupling a transport equation and an Hamilton-Jacobi equation.

##  2 Preliminary results

In this section, we report some of the results introduced in Bagagiolo-Pesenti
bagpes , and we present some new ways to approach the problem.

Suppose that a distributional evolution $t\mapsto{\cal M}(t)$ (the density of
the agents) is given. We denote by $V(\theta,t,w_{1},w_{2})$ the value
function (i.e. the infimum over the measurable controls of the cost functional
(2)), faced by an agent that at time $t$ is in the position
$(\theta,w_{1},w_{2})$ and starts its evolution. We drop the dependence of $V$
on $\cal M$. Assuming that the functions $t\mapsto F^{w_{1},w_{2}}({\cal
M}(t))$ are all continuous and bounded, then $V:B\times[0,T]\to\mathbb{R}$ can
be uniquely identified by the viscosity solutions of system of Hamilton-Jacobi
equations (one per each branch) with boundary conditions mutually exchanged
between the branches. On the switching points of the branch, the boundary
conditions are given by the value function itself evaluated on the point where
we switch on, in the new branch. This fact comes from the dynamic programming
principle and for a suitable interpretation of the optimal control problem on
a single branch as a suitable exit-time optimal control problem (see
Bagagiolo-Danieli bagdan , and the references therein for similar problems.).

Figure 2: Switching cross branches representation

Indeed, when the distribution $\cal M$ is given (i.e. when the fixed point
procedure described above is performed), on every branch (as in Figure 2) the
optimal control problem faced by any agent is a rather standard optimal
control given by some suitable combinations of reaching target (the station)
problem and exit time problems. In particular, on every branch
$B_{w_{1},w_{2}}$ we have the $(\theta,t)$-depending continuous value function
$V(\cdot,\cdot,w_{1},w_{2})$. The exit cost from $B_{1,1}$ at time $\tau$ is
given by $V(\theta_{1},\tau,0,1)$ in $\theta_{1}$ and by
$V(\theta_{2},\tau,1,0)$ in $\theta_{2}$. The exit costs from $B_{0,1}$ and
from $B_{1,0}$ are, respectively, $V(\theta_{2},\tau,0,0)$ and
$V(\theta_{1},\tau,0,0)$ (on both exit points of the branch). On the final
branch $B_{0,0}$ we only have the problem of reaching the station by the final
time $T$ paying as less as possible. Note that, when $\cal M$ is given, the
optimal control problem on $B_{0,0}$ can be solved (and in particular
$V(\cdot,\cdot,0,0)$ can be calculated) independently from the behavior on the
other branches. Solved the problem on $B_{0,0}$ we can then solve the problem
on $B_{0,1}$ and on $B_{1,0}$ and eventually on $B_{1,1}$.

Once the control problem is solved on every branch, then we tackle the optimal
transport problem which, when the initial distribution is given, returns the
evolution of the distribution of the agents in the interval $[0,T]$. We can
then perform the fixed point analysis.

More formally, for every branch $B_{w_{1},w_{2}}$, the optimal feedback
control $u^{*}$ is, at least in principle,

$u^{*}(\theta,t,w_{1},w_{2})=-V_{\theta}(\theta,t,w_{1},w_{2}),\quad\forall
t,\forall\theta.$ (4)

For every $(w_{1},w_{2})\in\\{0,1\\}\times\\{0,1\\}$ and for every
$t\in[0,T]$, we indicate by $m^{w_{1},w_{2}}(\cdot,t):[0,2\pi]\to[0,+\infty[$
the time dependent measure that represents the agents’ distribution on the
branch $B_{w_{1},w_{2}}$ at time $t$. We assume that the city network is
initially empty of tourist, that is $m^{w_{1},w_{2}}(\cdot,0)=0$ for
$(w_{1},w_{2})\in\\{0,1\\}\times\\{0,1\\}$. By conservation of mass principle,
$\cal M$ in point satisfies

$\displaystyle\int_{B}d{\cal M}(t)=$
$\displaystyle\int_{B_{1,1}}dm^{1,1}(t)+\int_{B_{1,0}}dm^{1,0}(t)+$
$\displaystyle+\int_{B_{0,1}}dm^{0,1}(t)+\int_{B_{0,0}}dm^{0,0}(t)=\int_{0}^{t}g(s)ds\quad
t\in[0,T].$

There a transport equations for the density $m^{w_{1},w_{2}}$ must hold, as an
example in case of $B_{1,0}$ we have:

$\displaystyle
m^{1,0}_{t}(\theta,t)-[V_{\theta}(\theta,t,1,0)m^{1,0}(\theta,t)]_{\theta}=0\
\mbox{in }B_{1,0}\times[0,T]$ $\displaystyle
m^{1,0}(\theta_{2},t)=m^{1,1}(\theta_{2},t),.$ (5)

The switching nature of the problem is reflected in the boundary conditions
like (5), which tie the distributions $m^{w_{1},w_{2}}$ of the different
branches $B_{w_{1},w_{2}}$.

As in bagpes , there are some suitable hypotheses: the initial distribution is
null (no-one is around the city at $t=0$); the congestion functions ${\cal
F}^{w_{1},w_{2}}$ do not explicitly depend on $\theta$ (all agents in the same
branch at the same instant equally suffer the same congestion pain). These
hypotheses lead to fact that the control choice each agent makes in one of the
point $(\theta_{S};1,1)$, $(\theta_{1},0,1)$, $(\theta_{2},1,0)$,
$(\theta_{1},0,0)$ and $(\theta_{2},0,0)$ do not change until the agent
remains in the same branch, and it is constant in time. Hence, only those
points are to be considered in the optimization procedure. In the sequel we
will call them: significant points.

Now here, instead of treating coupled Hamilton-Jacobi equations and transport
equation as it is common in mean field game theory), we write equivalent
conditions but in a different and in some sense more operative ways.

An agent in $(\theta_{1},0,0)$, at time $t\in[0,T]$, has two possibilities of
behavior (i.e. the optimal behavior is certainly one of the following two):
either to stay at $\theta_{1}$ without moving or to move and reach
$\theta_{S}$ by time $T$ (it is certainly not convenient to reach $\theta_{S}$
before the time $T$ and wait there; see the considerations made in bagpes ).
In the first case the optimal control is $u\equiv 0$, in the second case it is
$u\equiv\pm\frac{\theta_{S}-\theta_{1}}{T-t}$. Here $\theta_{s}-\theta_{1}$
stays for the length of the minimal path between $\theta_{S}$ and $\theta_{1}$
in the circular network, and the sign is coherently chosen with (1), in order
to run that path. Similarly when the agent is in $(\theta_{2},0,0)$.

Hence we get, given the cost functional (2),

$\begin{array}[]{ll}\displaystyle
V(\theta_{1},t,0,0)=\min\left\\{c_{S},\frac{1}{2}\frac{(\theta_{S}-\theta_{1})^{2}}{T-t}\right\\}+\int_{t}^{T}{\cal
F}^{0,0}\,ds\\\ \displaystyle
V(\theta_{2},t,0,0)=\min\left\\{c_{S},\frac{1}{2}\frac{(\theta_{S}-\theta_{2})^{2}}{T-t}\right\\}+\int_{t}^{T}{\cal
F}^{0,0}\,ds\end{array}$ (6)

where we do not display the argument of ${\cal F}^{0,0}$, being it fixed (the
distribution $\cal M$).

When the agent starts from $(\theta_{1},0,1)$ at the time $t$, it has three
choices: to stay at $\theta_{1}$ without moving; to reach $\theta_{S}$ at the
time $T$; to reach $\theta_{2}$ at some $\tau\in]t,T]$. Similarly as before,
in the first case it uses the control $u\equiv 0$; in the second case
$u\equiv\pm\frac{\theta_{S}-\theta_{1}}{T-t}$; and in the third case
$u\equiv\pm\frac{\theta_{2}-\theta_{1}}{\tau-t}$. Similarly for
$(\theta_{2},1,0)$. Hence:

$\displaystyle V(\theta_{1},t,0,1)=\min_{\tau\in[t,T]}$
$\displaystyle\left\\{c_{2}+c_{S}+\int_{t}^{T}{\cal
F}^{0,1}\,ds,c_{2}+\frac{1}{2}\frac{(\theta_{S}-\theta_{1})^{2}}{T-t}+\int_{t}^{T}{\cal
F}^{0,1}\,ds,\right.$
$\displaystyle\left.\frac{1}{2}\frac{(\theta_{2}-\theta_{1})^{2}}{\tau-t}+\int_{t}^{\tau}{\cal
F}^{0,1}\,ds+V(\theta_{2},\tau,0,0)\right\\}$ $\displaystyle\displaystyle
V(\theta_{2},t,1,0)=\min_{\tau\in[t,T]}$
$\displaystyle\left\\{c_{1}+c_{S}+\int_{t}^{T}{\cal
F}^{1,0}\,ds,c_{1}+\frac{1}{2}\frac{(\theta_{S}-\theta_{2})^{2}}{T-t}+\int_{t}^{T}{\cal
F}^{1,0}\,ds,\right.$
$\displaystyle\left.\frac{1}{2}\frac{(\theta_{1}-\theta_{2})^{2}}{\tau-t}+\int_{t}^{\tau}{\cal
F}^{1,0}\,ds+V(\theta_{1},\tau,0,0)\right\\}$ (7)

Finally, in $(\theta_{S},1,1)$ at the time $t$, the possibilities are: to stay
there without moving, to reach $\theta_{1}$ at a certain $\tau\in[t,T]$; and
to reach $\theta_{2}$ at a certain $\eta\in[t,T]$. Hence, similarly as before,

$\displaystyle\displaystyle V(\theta_{S},t,1,1)$
$\displaystyle=\min_{\tau,\eta\in[t,T]}\left\\{c_{1}+c_{2}+\int_{t}^{T}{\cal
F}^{1,1}\,ds,\frac{1}{2}\frac{(\theta_{1}-\theta_{S})^{2}}{\tau-t}+\int_{t}^{\tau}{\cal
F}^{1,1}\,ds+\right.$
$\displaystyle\left.+V(\theta_{1},\tau,0,1),\frac{1}{2}\frac{(\theta_{2}-\theta_{S})^{2}}{\eta-t}+\int_{t}^{\eta}{\cal
F}^{1,1}\,ds+V(\theta_{2},\eta,1,0)\right\\}$ (8)

Note that from (6)–(2), when $\cal M$ is fixed, performing the minimization
procedures we get an optimal behavior: to stay or to move towards some of the
other points. In doing that we also detect the possible arrival time ($\tau$
and $\eta$), and finally the corresponding optimal control $u$.

##  3 Existence of a mean field game equilibrium

In this section we give a proof of the existence of a mean field game
equilibrium.

By the considerations of the previous section, the only significant points are
$(\theta_{S},1,1)$, $(\theta_{1},0,1)$, $(\theta_{2},1,0)$, $(\theta_{1},0,0)$
and $(\theta_{2},0,0)$. Together with them, we have also to consider the flows
of arrivals on them (see Figure 2): $g$ the external arrival flow at the
station (it is a datum); $g_{01}$ the arrival flow in $(\theta_{1},0,1)$;
$g_{10}$ the arrival flow in $(\theta_{2},1,0)$; $g_{12}$ the arrival flow in
$(\theta_{2},0,0)$; $g_{21}$ the arrival flow in $(\theta_{1},0,0)$ (the last
four being part of the solution). Such flows must satisfy the conservation
conditions, for all $t\in[0,T]:$

$\begin{array}[]{ll}\displaystyle\int_{0}^{t}g(\tau)d\tau\geq\int_{0}^{t}g_{01}(\tau)d\tau+\int_{0}^{t}g_{10}(\tau)d\tau,\\\
\displaystyle\int_{0}^{t}g_{01}(\tau)d\tau\geq\int_{0}^{t}g_{12}(\tau)d\tau,\\\
\displaystyle\int_{0}^{t}g_{10}(\tau)d\tau\geq\int_{0}^{t}g_{21}(\tau)d\tau.\end{array}$
(9)

Denoting by $\rho^{w_{1},w_{2}}(t)$ the actual total mass of agents on the
branch $B_{w_{1},w_{2}}$, we then have

$\begin{array}[]{ll}\displaystyle\rho^{1,1}(t)=\int_{0}^{t}g(\tau)d\tau-\int_{0}^{t}g_{01}(\tau)d\tau-\int_{0}^{t}g_{10}(\tau)d\tau,\\\
\displaystyle\rho^{0,1}(t)=\int_{0}^{t}g_{01}(\tau)d\tau-\int_{0}^{t}g_{12}(\tau)d\tau,\\\
\displaystyle\rho^{1,0}(t)=\int_{0}^{t}g_{10}(\tau)d\tau-\int_{0}^{t}g_{21}(\tau)d\tau,\\\
\displaystyle\rho^{0,0}(t)=\int_{0}^{t}g_{12}(\tau)d\tau+\int_{0}^{t}g_{21}(\tau)d\tau.\end{array}$
(10)

Actually, using the notation of the previous section, it is

$\rho^{w_{1},w_{2}}(t)=\int_{B_{w_{1},w_{2}}}m^{w_{1},w_{2}}(\theta,t)d\theta.$

Form now on, we suppose that ${\cal F}^{w_{1},w_{2}}$ only depends, in a
continuous manner, from $\rho^{w_{1},w_{2}}$ and denote by $\rho$ the 4-uple
$(\rho^{1,1},\rho^{0,1},\rho^{1,0},\rho^{0,0})$. However, more general
situations may be taken into account. Hence, for what concerns our problem, to
give the distribution $\cal M$ is equivalent to give $\rho$. Accordingly, the
idea is to perform a fixed-point procedure on $\rho$ sent to a new
$\rho^{\prime}$ given by the transport law given by the minimization (6)–(2)
with $\rho$ inserted on it. To this end, we need some convexity and
compactness properties that we introduce in the following.

First of all note that the flow functions $g_{\cdot,\cdot}$ are bounded. Their
bounds come from the bound of the external arrival flow $g$ at $\theta_{S}$
(see also the Appendix for more details on such flows). Then functions
$\rho^{w_{1},w_{2}}$ are Lipschitz continuous, and moreover they are also
bounded by $K=\int_{0}^{T}g(s)ds$. Hereinafter, we denote by $L$ the Lipschitz
constant of any Lipschitz function we are going to consider. We then have, for
the quadruplet $\rho$,

$\rho\in\left\\{f:[0,T]\to[0,+\infty[\ \mbox{s.t. Lipschitz constant }\ L,\
\mbox{bounded by}\ K\right\\}^{4}:=X.$

Space $X$ is convex and compact for the uniform topology.

Next, obserserve that the function for which we need to exhibit the existence
of a fixed point should be a continuous function $\psi:X\to X$ which acts in
the following way: it takes $\rho\in X$, inserts $\rho$ in (6)–(2), obtains
the optimal controls $u$, constructs the flow-functions $g_{\cdot,\cdot}$ and
gives the new corresponding $\rho^{\prime}\in X$. Note that the flow-functions
$g_{01},g_{10},g_{12}$ and $g_{21}$ as well as the external flow $g$ are time
densities which, in order to run along the branches $B_{w_{1},w_{2}}$ with the
optimal control, at first must transform in spatial densities and then become
again time densities on the switching points. Hence, the function $\psi$ need
to built the spatial and temporal components of the mentioned flow-functions
in order to obtain the new $\rho^{\prime}\in X$ (we report the computation of
both spatial and temporal components of the above flow-functions in the
Appendix.)

As regards the continuity of $\psi$ the main difficulty is that the optimal
feedback controls may be not unique: for any fixed $t$, the minimization
procedure in (6)–(2) may return more than one minimizers, when more than one
of the arguments of the considered $\min\\{.,.\\}$ functions reaches the
minimum value. In particular, this situation may happen for several $t$, even
accumulating or in a whole interval, in (7)–(2). Differently, in (6) it may
occur at most in a single instant and in this case is not a problem.

When such a multiplicity situation occurs, then the agents naturally split in
more communities using one of the optimal behaviors. In some sense they use
mixed-strategies, that is they convexify. The function $\psi$ described above,
is defined assuming that the agents do not split: all agents, in the same
position and at the same instant makes the same choice for their behavior
(this is indeed a characteristic feature of the mean field games model: the
agents are homogeneous and indistinguishable). But this is coherent only with
the uniqueness of the optimal feedback. When, instead, we have multiplicity,
this is not more true and, when such a multiplicity accumulates in time, the
continuity of the function $\psi$ is not evident. In particular there are more
than one way to construct it: just taking different optimal behaviors.

In Bagagiolo-Pesenti bagpes , the above difficulty is overcome by a-priori
assuming that there exists a finite number $N$ that bounds the maximum number
of times at which those multiplicities appear, independently on $\rho$. This
assumption guarantees a sort of continuity of the function $\psi$ by
considering a multi-function whose image is the closed convex hull of all
possible outputs of the function $\psi$. Specifically, the assumption on the
uniform bound $N$, guarantees that such a closed convex hull is given by a set
of a a-priori bounded number of extremal points, and hence getting compact
image and closed graph.

Here we let such an assumption drop and we prove the continuity and convexity
of $\psi$ by defining a sequence of continuous convex functions approximating
of $\psi$, based on the assumption that an agent can make sub-optimal
decisions as described in the following.

Take $\varepsilon>0$ and, whenever you make an optimal choice for all agents
in one of the minimization problems (6)–(2), then maintain the same choice,
i.e., consider the same argument within the $\min\\{.,.\\}$ function as
minimizing one until, possibly, its difference with a new minimizing one is
greater than $\varepsilon$. In this latter case, select the new minimizer as
the new choice. Calling $\psi_{\varepsilon}$ such a new correspondence we now
describe it in a precise way.

In (2), let us call (2)-1, (2)-2, (2)-3 respectively the first, the second and
the third term inside the parenthesis in the right-hand side. Note that, when
(2)-1 is the minimizer, then it means that to stay at $\theta_{S}$ without
moving till the final time $T$ is an optimal strategy. When the minimum is
realized by (2)-2 (for some $\tau$), then to move and reach $\theta_{1}$ at
time $\tau$ is an optimal strategy. When the minimum is given by (2)-3 (for
some $\eta$), then to move and to reach $\theta_{2}$ at time $\eta$ is an
optimal strategy. In the sequel, when we say, for example, “at the time $t$
the agents make the choice (2)-2”, it means that they also choose the
corresponding minimizing $\tau$ and move towards $\theta_{1}$ with constant
velocity $|\theta_{1}-\theta_{S}|/(\tau-t)$. But it does not mean that it is
necessary the best choice at all.

At $(\theta_{S},1,1)$ at the time $t_{0}=0$, all the present agents make the
same optimal choice, among the optimal ones, for example (2)-2. Then, all
other agents, at $(\theta_{S},1,1)$ at every subsequent times $t$, continue to
make the same choice (2)-2, as long as, at $t$, (2)-2 realizes the minimum up
to an error lower than $\varepsilon$ (i.e. they continue to choose (2)-2 even
if it is not more optimal, but anyway $\varepsilon$-optimal). Let
$t_{1}\in[0,T]$ be the possible first instant at which the error is equal to
$\varepsilon$. Then the agents at $(\theta_{S},1,1)$ at $t_{1}$ make a new
optimal choice (certainly different from (2)-2) among the optimal ones, for
example (2)-3. Then, all other agents, at $(\theta_{S},1,1)$ at every
subsequent times $t$, continue to make the same choice (2)-3, as long as, at
$t$, (2)-3 realizes the minimum up to an error lower than $\varepsilon$. We
proceed in this way, taking a possible new instant $t_{2}$, and so on.

We similarly argue for the other points $(\theta_{2},1,0)$,
$(\theta_{1},0,1)$, $(\theta_{2},0,0)$, $(\theta_{1},0,0)$ and for the whole
time interval $[0,T]$. Note that the arguing above is not influenced by the
actual presence of the agents in the points.

Of course, starting from $t_{0}=0$, we may have more than one of such possible
constructions. For example, at $(\theta_{S},1,1)$ we may have more than one
possible optimal choices at the time $t_{0}=0$. Different choices will
generate different first instants $t_{1}$ for which other optimal choices will
be considered, which will generate different second instants $t_{2}$ and so
on. However, the number of all such possible instants is finite. Indeed, the
functions involved in the minimization (6)–(2) can be considered as Lipschitz
functions of their time-entries: $t$, $\tau$, $\eta$. Actually, the terms of
the form, for example, $(\theta_{1}-\theta_{2})^{2}/(\tau-t)$ are not
Lipschitz in $t,\tau\in[0,T]$, but they play a role in the minimization
certainly until they are not grater than $c_{s}+c_{1}+c_{2}$, and until it
happens, they are uniformly Lipschitz. Hence, at every significant points as
above, the distance of two subsequent instants in anyone of the sequences
constructed as before is not less than $\varepsilon/L$. This means that, in
the time interval $[0,T]$, we have a uniform a-priori bound $N_{\varepsilon}$
for the number of all those possible instants. Hereinafter $N_{\varepsilon}$
will denote any possible bound of the number of some specific quantities.

The above arguments allow us to state that for each one of the significant
points as above we can divide the interval $[0,T]$ in a finite number of sub-
intervals, in the interior of which the number and the type of the
$\varepsilon$-optimal choices does not change. Such choices are indeed
$\varepsilon$-optimal because they realize the minimum in the corresponding
formula (6)–(2) up to an $\varepsilon$-error, and those minima are exactly the
optima.

Now we describe as $\psi_{\varepsilon}$ acts. Take
$\rho=(\rho^{1,1},\rho^{0,1},\rho^{1,0},\rho^{0,0})\in X$ and insert it in
(6)–(2); for each significant point construct all possible sequences of
$\varepsilon$-optimal behaviors choosing one of them inside anyone of the time
sub-intervals (the number of all such sequences, as well as the number of
their terms is bounded by $N_{\varepsilon}$). For each one of those sequences,
starting by the given external arrival flow $g$, construct the flow functions
$g_{\cdot,\cdot}$ (see the Appendix) and then the new mass concentration
$\rho^{\prime}=(\rho^{\prime 1,1},\rho^{\prime 0,1};\rho^{\prime
0,1},\rho^{\prime 0,0})\in X$, using (10). Define
$\psi_{\varepsilon}(\rho)\subseteq X$ as the set of all possible mass
concentrations $\rho^{\prime}$ constructed in this way. Note that
$\psi_{\varepsilon}(\rho)$ is a finite set with no more than $N_{\varepsilon}$
elements, independently on $\rho$.

Since $\psi_{\varepsilon}$ is a multifunction and we want to apply the fixed-
point Kakutani Theorem, we need compact and convex images and closed graph.
The set $\psi_{\varepsilon}(\rho)$ being finite is certainly not convex, hence
we have to convexify. We define as $\tilde{\psi}_{\varepsilon}(\rho)\subseteq
X$ the set of all quadruplets constructed in the following way. At
$(\theta_{S},1,1)$, for each one of the sub-intervals of $[0,T]$ suppose that
the agents (arriving with the given flow $g$) split in some fractions each one
of them following one of the possible $\varepsilon$-optimal choices in the
time sub-intervals. This gives a mass concentration $\tilde{\rho}^{1,1}$ on
the branch $B_{1,1}$, and, considering all the admissible splitting, we get
the first components of all elements of $\tilde{\psi}_{\varepsilon}(\rho)$.
For any of the possible splitting we also get the corresponding flow functions
$g_{01},g_{10}$, and then, as above, we consider all admissible splitting in
the points $(\theta_{1},0,1)$ and $(\theta_{2},1,0)$ getting the mass
concentrations in the corresponding branches. This gives the second and the
third components of the elements of $\tilde{\psi}_{\varepsilon}(\rho)$. In a
similar way we construct the fourth components. The set
$\tilde{\psi}_{\varepsilon}(\rho)$ is then a convex subset of $X$ (any split
is a convex combination of masses). It is a sort of convex hull of
$\psi_{\varepsilon}(\rho)$, suitably performed branch by branch. The fact that
$\psi_{\varepsilon}(\rho)$ is a finite set, also gives the closedness (and
hence the compactness, being $X$ compact) of
$\tilde{\psi}_{\varepsilon}(\rho)$ in $X$ (for the uniform topology).
Similarly, the multifunction $\rho\mapsto\tilde{\psi}_{\varepsilon}(\rho)$ has
closed graph. Here, we sketch such a proof.

We have to prove that whenever $\rho_{n}$ converges to $\rho$ in $X$, and
$\xi^{n}\in\tilde{\psi}_{\varepsilon}(\rho^{n})$ converges to $\xi$ in $X$,
then $\xi\in\tilde{\psi}_{\varepsilon}(\rho)$. By definition, the first
component of the quadruplet $\xi^{n}$ is given by a suitable split of the
external arrival flow $g$ in every sub-intervals of $[0,T]$ and hence by the
transport with the corresponding choices (2)-1, (2)-2, (2)-3. Since,
independently on $\rho^{n}$, the number of sub-intervals is bounded by
$N_{\varepsilon}$, and the number of choices is always not more than three (in
other words, the number of the element of $\psi_{\varepsilon}(\rho^{n})$ is
uniformly bounded with respect to $\rho$), we may suppose that, at least for a
subsequence, the first component of all $\xi^{n}$ is constructed by using the
same number of sub-intervals with exactly the same type of choices for every
sub-interval. Moreover, the sub-interval depending triple,
$(\lambda^{n}_{1},\lambda^{n}_{2},\lambda^{n}_{3})$ representing the sequences
of split-fractions of the mass of agents that choose one of the admissible
choices (if in the sub-interval a choice is not admissible, then the
corresponding $\lambda$ is zero) converge, sub-interval by sub-interval, to a
limit fraction $(\lambda_{1},\lambda_{2},\lambda_{3})$. Due to the convergence
of $\rho^{n}$ to $\rho$, using (2) we have that the first sub-interval of
every step $n$ converges to a limit sub-interval, which is coherent with (2)
and $\rho$. Similarly, for the second sub-intervals and so on. Hence we have a
candidate limit partition of $[0,T]$ in sub-intervals, with the corresponding
choices and the corresponding split-fractions. Since $\xi^{n}$ uniformly
converges to $\xi$, it is easy to see that the first component of $\xi$ is
generated by those sub-intervals, choices and split-fractions. Moreover, sub-
interval by sub-interval, the flow functions to the branches $B_{0,1}$ and
$B_{1,0}$ converges to some some flow functions which, reasoning as above in
the new branches gives the construction (sub-intervals, choices, split-
fraction) of the second and third component of $\xi$. Similarly for the fourth
component. Hence $\xi\in\tilde{\psi}_{\varepsilon}(\rho)$.

By the Kakutani fixed point theorem, for every $\varepsilon>0$ there exists a
fixed point $\rho_{\varepsilon}\in X$ for $\tilde{\psi}_{\varepsilon}$, that
is $\rho_{\varepsilon}\in\tilde{\psi}_{\varepsilon}(\rho_{\varepsilon})$. This
means that, if all agents suppose that the realized congestion is given by
$\rho_{\varepsilon}$ than there is a suitable split choice of behaviors such
that $\rho_{\varepsilon}$ is actually realized and moreover, by construction,
every agent, when moving by one of the admissible choices (split fractions not
null), is subject to a cost (2) whose difference with the optimal one is
infinitesimal as $\varepsilon$. In some sense $\rho_{\varepsilon}$ is a
$\varepsilon$\- mean field games equilibrium.

###### Definition 1

A mean field game equilibrium for our problem is a pair $(\rho,\lambda)$ where
$\rho\in X$ and $\lambda\in L^{\infty}(0,T)^{13}$ is a string of $13$ split-
fractions functions, that is

$\begin{array}[]{ll}\displaystyle\lambda=\left(\lambda^{(\theta_{S},1,1)}_{1},\lambda^{(\theta_{S},1,1)}_{2},\lambda^{(\theta_{S},1,1)}_{3},\right.\\\
\displaystyle\lambda^{(\theta_{1},0,1)}_{1},\lambda^{(\theta_{1},0,1)}_{2},\lambda^{(\theta_{1},0,1)}_{3},\lambda^{(\theta_{2},1,0)}_{1},\lambda^{(\theta_{2},1,0)}_{2},\lambda^{(\theta_{2},1,0)}_{3},\\\
\displaystyle\left.\lambda^{(\theta_{1},0,0)}_{1},\lambda^{(\theta_{1},0,0)}_{2},\lambda^{(\theta_{2},0,0)}_{1},\lambda^{(\theta_{2},0,0)}_{2}\right)\end{array}$
(11)

which means, for example, that in $(\theta_{S},1,1)$ we consider an arrival
flow of density $\lambda^{(\theta_{S},1,1)}_{1}g$ of agents which will make
the choice (2)-1, an arrival flow of density $\lambda^{(\theta_{S},1,1)}_{2}g$
of agents which will make the choice (2)-2, and so on. This will generate
flows functions on the branches $B_{0,1}$ and $B_{1,0}$ which will be
correspondingly split by the fractions in the second line of (11), and so on.

The pair $(\rho,\lambda)$ is a mean field game equilibrium if the behavior
explained above generates exactly the mass density $\rho$ (by (10)) and if,
putting $\rho$ inside the cost (2), that behavior makes each agent optimize.
This means that for almost every instants $t$, whenever one of the split
fractions is not null, then the corresponding choice is an optimal one.

###### Remark 1

Similarly to the definition of $\lambda$ for the mean field games equilibrium
(11), we have the split fractions functions $\lambda_{\varepsilon}$ for the
$\varepsilon$-problem. In that kind of problem, such functions are, by
construction, constant in the sub-intervals. Denoting by
$\overline{\lambda}_{\varepsilon}$ the string of $13$ of split functions we
used in the construction of the element of $\psi_{\varepsilon}(\rho)$ as well
as of $\tilde{\psi}_{\varepsilon}(\rho)$, the corresponding notion of split
functions $\lambda_{\varepsilon}$ as in (11) is given, in some sense, by the
product of elements of $\overline{\lambda}_{\varepsilon}$. For example

$\begin{array}[]{ll}\displaystyle(\lambda_{\varepsilon})_{2}^{(\theta_{S},1,1)}=(\overline{\lambda}_{\varepsilon})_{2}^{(\theta_{S},1,1)}\\\
\displaystyle(\lambda_{\varepsilon})_{2}^{(\theta_{1},0,1)}=(\overline{\lambda}_{\varepsilon})_{2}^{(\theta_{S},1,1)}(\overline{\lambda}_{\varepsilon})_{2}^{(\theta_{1},0,1)}\\\
\displaystyle\mbox{(the product between the arriving fraction and the
departing fraction)}\\\
\displaystyle(\lambda_{\varepsilon})_{2}^{(\theta_{1},0,0)}=(\overline{\lambda}_{\varepsilon})_{2}^{(\theta_{S},1,1)}(\overline{\lambda}_{\varepsilon})_{2}^{(\theta_{1},0,1)}(\overline{\lambda}_{\varepsilon})_{2}^{(\theta_{1},0,0)}\\\
\displaystyle\mbox{(here the arriving fraction is the product of the first two
factors)}\end{array}$ (12)

All the fractions functions must applied to the global flow $g$ transported,
branch by branch, by the chosen control. We refer to the Appendix for more
details in the subject of transport on networks.

###### Theorem 3.1

There exists a mean field game equilibrium

Proof. For every $\varepsilon>0$ we have a fixed point $\rho_{\varepsilon}$
for $\tilde{\psi}_{\varepsilon}$, and we have the corresponding split-
fractions function $\lambda_{\varepsilon}\in L^{\infty}(0,T)^{13}$ which, for
every $\varepsilon$ is constant in each sub-intervals (note that the numbers
of sub-intervals is not in general bounded uniformly in $\varepsilon$). Let
$\lambda\in L^{\infty}(0,T)^{13}$ be the weak-star limit of
$\lambda_{\varepsilon}$ as $\varepsilon\to 0$ (i.e.
$\int_{0}^{t}\lambda_{\varepsilon}f\to\int_{0}^{T}\lambda f$ for all
integrable function $f$) and $\rho\in X$ be the uniform limit of
$\rho_{\varepsilon}$ (at least for a subsequence). By the construction of
$\lambda_{\varepsilon}$, by the weak-star convergence of
$\lambda_{\varepsilon}$ to $\lambda$, and by the uniform convergence of
$\rho_{\varepsilon}$ to $\rho$, we get that $\lambda$ is a split-fractions
function generating the mass density $\rho$ (it generates flow functions
$g_{\cdot,\cdot}$ which by (10) generate $\rho$). Finally, at every step
$\varepsilon$, if the agents exactly perform the behavior given by the split-
fractions $\lambda_{\varepsilon}$, then they get an $\varepsilon$-optimal
cost. We then conclude that the pair $(\rho,\lambda)$ is a mean field game
equilibrium. Indeed, for almost every instants $t$ such that one of the
components of $\lambda$ is not null, by the weak-star convergence there is, at
least for a sub-sequence, a sequence of instants $t_{\varepsilon}$ converging
to $t$ such that the corresponding component of
$\lambda_{\varepsilon}(t_{\varepsilon})$ is not null. This means that the
corresponding choice is $\varepsilon$-optimal with respect to
$\rho_{\varepsilon}$ and we conclude. ∎

###### Remark 2

The introduction of the $\varepsilon$-problem, is mostly due to an analytical
purpose: the possibility of the construction of an $\varepsilon$-mean field
games equilibrium
$\rho_{\varepsilon}\in\tilde{\psi}_{\varepsilon}(\rho_{\varepsilon})$ by the
fact that the number of sub-intervals is finite, and then the passage to the
limit.

However, from an applicative modeling point of view, the $\varepsilon$-problem
is also interesting by itself. Indeed, it takes account of a typical human
decision phenomenon: herd behavior (see, e.g., Banerjee ), where an agent
looks at and is influenced by the decisions made by other agents in taking its
own decision, even if it is not the optimal choice for him. In these
situations a second thought occurs when discrepancy from the optimal choice
becomes too large (more than $\varepsilon$) and than it is evident that it is
better to change choice.

##  4 The optimization problem

In this section we introduce a possible optimization problem faced by a local
authority, hereinafter referred as to controller, that intends to manage the
flow of tourists.

We restrict our analysis to congestion cost functions of the form

${\cal
F}^{w_{1},w_{2}}(\rho)=\alpha_{w_{1},w_{2}}\rho^{w_{1},w_{2}}(s)+\beta_{w_{1},w_{2}}$
(13)

with $\rho=(\rho^{1,1},\rho^{0,1},\rho^{1,0},\rho^{0,0})\in X$. In (13), the
coefficients $(\alpha_{w_{1},w_{2}},\beta_{w_{1},w_{2}})$ are at disposal of
an external controller whose aim is to force the realized mean field game
string $\rho$ to be as close as possible (in the uniform topology) to a
reference string $\overline{\rho}\in X$, i.e. to minimize:

$\max_{w_{1},w_{2}\in\\{0,1\\}}\left\\{\max_{t\in[0,T]}|\rho^{w_{1},w_{2}}(t)-\overline{\rho}^{w_{1},w_{2}}(t)|\right\\}=\|\rho-\overline{\rho}\|_{X}$
(14)

Formally, let us denote by $\chi_{\alpha,\beta}$ the set of all possible (mass
component) mean field games equilibrium corresponding to the choice of the
parameters $\alpha=(\alpha_{1,1},\alpha_{0,1},\alpha_{1,0},\alpha_{0,0})$ and
$\beta=(\beta_{1,1},\beta_{0,1},\beta_{1,0},\beta_{0,0})$, and let us suppose
that the parameters are restricted to belong to a compact set
$K\subset\mathbb{R}^{4}\times\mathbb{R}^{4}$. Then, the controller faces the
optimization problem given by

$\inf_{(\alpha,\beta)\in K}{\cal E}(\alpha,\beta)=\inf_{(\alpha,\beta)\in
K}\left(\inf_{\rho\in\chi_{\alpha,\beta}}\|\rho-\overline{\rho}\|_{X}\right)$
(15)

###### Theorem 4.1

There exists an optimal pair $(\alpha,\beta)\in K$.

Proof. Let $(\alpha^{n},\beta^{n})\in K$ be a minimizing sequence for $\cal
E$, and for every $n$ let $\rho^{n}\in\chi_{\alpha_{n},\beta_{n}}$ realize the
infimum in (15) up to $1/n$. By compactness we may suppose that
$(\alpha^{n},\beta^{n})$ converges to $(\tilde{\alpha},\tilde{\beta})\in K$,
and that $\rho^{n}$ uniformly converges to $\tilde{\rho}\in X$. To prove the
statement, we only need to prove that
$\tilde{\rho}\in\chi_{\tilde{\alpha},\tilde{\beta}}$.

1) $\tilde{\rho}\in\chi_{\tilde{\alpha},\tilde{\beta}}$. This means to prove
that $\tilde{\rho}$ is a mean field game equilibrium for the problem. Let
$\lambda^{n}$ be the split-fraction functions associated to the mean field
game equilibrium $\rho^{n}$, and let $\tilde{\lambda}$ be a weak-star limit of
it. Then, due to the considered convergences, arguing as in the proof before,
we get that the pair $(\tilde{\rho},\tilde{\lambda})$ is a mean field
equilibrium for $(\tilde{\alpha},\tilde{\beta})$

∎

###### Remark 3

To consider in (15) an optimization problem of the form

$\inf_{(\alpha,\beta)\in
K}\left(\sup_{\rho\in\chi_{\alpha,\beta}}\|\rho-\overline{\rho}\|_{X}\right)$

would be certainly more interesting and robust from the point of view of an
external controller. However, the existence of an optimal pair, due to the
possibilities of more than one mean field game equilibrium is not evident.
Certainly there exists a pair $(\alpha,\beta)$ and $\rho\in{\cal
X}_{\alpha,\beta}$ such that $\|\rho-\overline{\rho}\|_{X}$ is equal to that
infimum, but this does not guarantee that the pair $(\alpha,\beta)$ is
optimal.

If we would be able, under some suitable hypotheses, to guarantee the
uniqueness of the mean field games equilibrium, then such a problem will be
bypassed. Maybe, stronger hypotheses on the costs $\cal F$ could be useful. We
leave such a question to future investigations.

So far we have assumed that the coefficients
$(\alpha_{w_{1},w_{2}},\beta_{w_{1},w_{2}})$ of congestion cost functions
${\cal F}^{w_{1},w_{2}}$ are not time-depending. This hypothesis plays an
essential role in the analysis of the problem, in particular on the nature of
the optimal controls. Here we little bit amend such a hypothesis by assuming
that $(\alpha_{w_{1},w_{2}},\beta_{w_{1},w_{2}})$ may be piecewise constant
functions and that the control is implemented at the significant points of
each branch of the network, e.g., through gates.

We consider a finite sequence of fixed instant
$t_{0}=0<t_{1}<t_{2}<...<t_{N}=T$, and for every $i=0,..N-1$, the coefficients
$(\alpha^{i},\beta^{i})\in K$ at disposal of the external controller in the
interval $[t_{i},t_{i+1}[$. The congestion cost paid by an agent at time $s$
becomes

$\tilde{\cal
F}^{w_{1},w_{2}}(s)=\alpha_{w_{1},w_{2}}(\tau(s))\rho^{w_{1},w_{2}}(s)+\beta_{w_{1},w_{2}}(\tau(s))$

where $\alpha_{w_{1},w_{2}}(\tau)=\alpha_{w_{1},w_{2}}^{i}$, respectively
$\beta_{w_{1},w_{2}}(\tau)=\beta_{w_{1},w_{2}}^{i}$, for
$t_{i}\leq\tau<t_{i+1}$, $i=0,...,N-1$; and $\tau(s)$ is the last switching
instant not greater than $s$ along the agent trajectory, i.e., the instant at
which the agent state entered $B_{w_{1},w_{2}}$.

In this situation the total cost payed by an agent is the usual

$\begin{array}[]{ll}J(t,u)=\int_{t}^{T}\left(\frac{u(s)^{2}}{2}+\tilde{\cal
F}^{w_{1}(s),w_{2}(s)}(s)\right)ds+\\\
c_{1}w_{1}(T)+c_{2}w_{2}(T)+c_{S}\xi_{\theta=\theta_{S}}(T)\end{array}$ (16)

but we have an explicit dependence on the time of the congestion cost through
$\alpha_{w_{1},w_{2}}(.)$ and $\beta_{w_{1},w_{2}}(.)$. Nevertheless, such a
dependence is compatible with the structure of the choices of the agents in
our model. Indeed, by the cost (16), the agents that arrive at
$(\theta_{S},1,1)$ at the beginning feel the actual value of the parameters
$\alpha^{1,1}$ and $\beta^{1,1}$, and act as they are constant until the
possibly exit from the branch $B_{1,1}$. Suppose that they switch on the
branch $B_{0,1}$. Then, at that moment, they feel the actual value of the
parameters $\alpha^{0,1}$ and $\beta^{0,1}$ and, again, act as they are
constant until the possibly exit from the branch, and so on. It is here worth
recalling, that we assume that the controller applies its controls at the
beginning of the network branches and that this latter model as the previous
one implies that the agents make their decisions each time they enter a new
branch, then a possible change in the parameters $(\alpha^{i},\beta^{i})$ is
not perceived by the agents that are “on the road” along a branch.
Consequently the considered time-dependence of the cost does not change the
reasoning we did before about properties of the optimal controls. In
particular, formulas (6)–(2) remain valid, and the existence of a mean field
games equilibrium is guaranteed. Moreover a similar optimization problems as
above can be successfully performed.

##  5 Conclusions

In this paper we have introduced a mean field game model as a possible way of
representing the flows of tourists in the alleys of an heritage city. Then, we
used this model as a basis of the optimization problem that aims at managing
the tourist flows according to some targets posed by the local authorities.

Many further refinements of both the model and the optimization problem
presented are possible from both a theoretical and an applicative perspective.
As an example, the assumption of the existence of surveillance cameras, to
count people entering and leaving areas of interest, may suggest the
description of a new type of information in the model and justify the
definition of a feedback control policy on the coefficients
$(\alpha_{w_{1},w_{2}},\beta_{w_{1},w_{2}})$ at disposal of an external
controller. In addition different kinds of local authority management measures
can be modeled and more than two places of interests may be considered. Note
that in the latter case the number of alternative paths would increase
exponentially. However the direct experience of the authors in Venice suggests
that the vast majority of the tourists is interested in very few attractions
in each city.

###### Acknowledgements.

This research was partially funded by a INdAM-GNAMPA project 2017

## References

  * (1) Achdou, Y., Camilli, F., Capuzzo Dolcetta, I.: Mean field games: numerical methods for the planning problem, SIAM J. Control Optim. 50(1), 77-109 (2012)
  * (2) Ambrosio, L., Gigli, N., Savarè, G.: Gradient Flows in Metric Spaces and in the Space of Probability Measures, Second Edition, Lectures in MAthematics, ETH Zürich, Birkhäuser Verlag, Basel, 2008.
  * (3) Bagagiolo, F., Danieli, K.:Infinite horizon optimal control problems with multiple thermostatic hybrid dynamics, Nonlinear Analysis; Hybrid Systems, 6 (2011), 824–838.
  * (4) Bagagiolo, F., Pesenti, R.: Non-memoryless Pedestrian Flow in a Crowded Environment with Target Sets, in: Apaloo J., Viscolani B. (eds) Advances in Dynamic and Mean Field Games. ISDG 2016. Annals of the International Society of Dynamic Games, vol 15. Birkhäuser, Cham, (2017)
  * (5) Bagagiolo, F., Bauso, D., Maggistro, R., Zoppello, M.: Game theoretic decentralized feedback controls in Markov jump processes, J Optim. Theory Appl., 173, 704-726 (2017)
  * (6) Banerjee, A. V.: A Simple Model of Herd Behavior, Q J Econ., 107(3), 797-817, (1992)
  * (7) “Brian” Park, B., Yun, I., Ahn, K.: Stochastic optimization for sustainable traffic signal control, International journal of sustainable transportation, 3(4), 263-284 (2009)
  * (8) Camilli, F., Carlini, E., Marchi, C.: A model problem for mean field games on networks, Discrete Contin. Dyn. Syst. 35, 4173-4192 (2015)
  * (9) Camilli, F., De Maio R., Tosin, A.: Transport of measures on networks, Netw. Heterog. Media 12, 191-215 (2017).
  * (10) Cardaliaguet, P.: Notes on mean field games. Unpublished notes.
https://www.ceremade.dauphine.fr/cardalia/index.html

  * (11) Carlier, G., Jimenez, C., Santambrogio, F.: Optimal transportation with traffic congestion and Wardrop equilibria, SIAM J. Control Optim., 47, 1330-1350, (2008)
  * (12) Carlier, G., Santambrogio, F.: A continuous theory of traffic congestion and Wardrop equilibria, J. Math Sci., 181 (6), 792-804, (2012)
  * (13) Cole, R., Dodis, Y., Roughgarden, T.: How much can taxes help selfish routing?, J. Comput. Syst. Sci., 72, 444-467 (2006)
  * (14) Como, G., Savla, K., Acemoglu, D., Dahleh, M. A., Frazzoli, E.: Stability analysis of transportation networks with multiscale driver decisions, SIAM J. Control Optim., 51(1), 230-252 (2013)
  * (15) Cristiani, E., Priuli, F., Tosin, A.: Modeling rationality to control self-organization of crowds: an environmental approach, SIAM J. Appl. Math. 75, 605-629 (2015)
  * (16) Dial, R. B.: Network-optimized road pricing: Part I: A parable and a model, Oper. Res., 47, 54-64 (1999)
  * (17) Gomes, G., Horowitz, R.: Optimal freeway ramp metering using the asymmetric cell transmission model, Transport Res C: Emer, 14(4), 244-262 (2006)
  * (18) Hegyi, A., De Schutter, B., Hellendoorn, H.: Model predictive control for optimal coordination of ramp metering and variable speed limits, Transport Res C: Emer 13(3), 185-209 (2005)
  * (19) Huang, M., Caines, P., Malhame, R.: Large population dynamic games; closed-loop McKean- Vlasov systems and the Nash certainly equivalence principle. Commun. Inf. Sys. 6, 221-252 (2006)
  * (20) Khattak, A., Polydoropoulou, A., Ben-Akiva, M.: Modeling revealed and stated pretrip travel response to advanced traveler information systems, Transport. Res. Rec., 1537, 46-54 (1996)
  * (21) Lachapelle, A., Salomon, J., Turinici, G.: Computation of mean-field equilibria in economics. Math. Models Methods Appl. Sci. 20(4), 1-22 (2010)
  * (22) Lasry, J.M., Lions, P.L.: Juex à champ moyen ii. Horizon fini et controle optimal. C. R. Math. 343, 679-684 (2006)
  * (23) Morrison, S.: A survey of road pricing, Transp. Res. A: Gen., 20(2), 87-97, (1986).
  * (24) Smith, M. J.: The marginal cost taxation of a transportation network, Transp. Res. B: Methodol., 13(3), 237-242 (1979)
  * (25) Srinivasan, K., Mahmassani, H.: Modeling inertia and compliance mechanisms in route choice behavior under real-time information, Transport. Res. Rec., 1725, 45-53 (2000)
  * (26) Wardrop, J. G.: Some theoretical aspects of road traffic research, Proc. Inst. Civ. Eng., 1(3), 325-362, (1952)

##  Appendix A Appendix: on some flows and transports on networks

In this appendix, we discuss the characteristics of the flow functions $g$ and
$g_{.,.}$ both from a spatial and a temporal perspective.

Le us consider what happens in our model on the branch $B_{1,1}$. At the point
$\theta_{S}$ there is an entering flow of agents, represented by the function
of time $g$. This is a “density with respect to time”, that is, from a
dimensional point of view, something like ”number of agents” divided by time.
When instead we calculate the actual total mass $\rho^{1,1}(t)$ on the branch,
we have to integrate on space (over the whole branch) a “density with respect
to space ”, that is something as ”number of agents” divided by space.
Moreover, in the two switching points $\theta_{1},\theta_{2}$ of the branch,
we need to recover again a “density with respect to time” which has to
represent the entering flow in the other branches. Actually, in our model,
being the congestion cost depending only on the total mass on the branch, we
do not really need to pass through the determination of density with respect
to space, but the computation of the exiting flows is enough, see (10).
However, in a more general setting we need both computations. In the sequel,
we call them the spatial and temporal component of the flow. Obviously, this
argument is linked to what is called, for example, disintegration of a measure
(see for instance Ambrosio-Gigli-Savaré ambgigsav  and Camilli-De Maio-Tosin
CDeMTos ).

Still considering the branch $B_{1,1}$ as example, we recall that, in our
model, the agents arriving in $\theta_{S}$ at time $t$ choose their behavior
at exactly that time and they do not change it until they, possibly, leave the
branch.

Now we sketch how, in our situation, to split the flow in the two components.
If an agent at time $t$ is in the position $\theta$ and it is going towards
$\theta_{1}$, then there is a time $\tau<t$ at which it has arrived at
$\theta_{S}$, and there is a time $\overline{\tau}>t$ at which it will arrive
at $\theta_{1}$. And it is using a constant velocity $u(\tau)$ of absolute
value $|\theta_{1}-\theta_{S}|/(\overline{\tau}-\tau)$. In particular here
$\overline{\tau}$ is calculated as the minimizing time for the second term in
(2) with $t$ replaced by $\tau$. Let us indicate by $\Lambda(\theta,t)$ the
function that gives the entering time $\tau$. Note that, when the mass
concentration is given (i.e. when performing the fixed point procedure) the
function $\Lambda$ is evaluable using (6)–(2) which give all possible constant
velocities $u(\tau)$.

In the present heuristic arguing we suppose that $\Lambda$ and $u$ are
differentiable (actually, in our setting they are Lipschitz and hence at least
almost everywhere differentiable). Assuming for simplicity $\theta_{S}=0$,
$\theta>0$, and denote $\lambda(\tau)=1/u(\tau)$, we have the following
relation

$\Lambda(\theta,t)+\lambda(\Lambda(\theta,t))\theta=t,$ (17)

Deriving with respect to $t$ and $\theta$ we have the relations

$\Lambda_{t}(\theta,t)=\frac{1}{1+\lambda^{\prime}(\Lambda(\theta,t))\theta}=-\frac{\Lambda_{x}(\theta,t)}{\lambda(\Lambda(\theta,t))}$

Now, supposing that all agents arriving in $\theta_{S}$ with flow $g$ move
towards $\theta_{1}$, then the density in time $g$ at $\theta_{S}$ is spread
on the spatial density component $-g(\Lambda(\theta,t))\Lambda_{x}(\theta,t)$
and on the temporal density component
$g(\Lambda(\theta,t))\Lambda_{t}(\theta,t)$ (which act only if in the point
$\theta$ at time $t$ actually there are agents, some agents are already
arrived, otherwise they are zero). Such relations may be verified by standard
mass balance/conservation arguments. In particular, at the point $\theta_{1}$
the arriving flow (density in time) which will be the entering flow in the new
branch $B_{0,1}$ is $t\mapsto
g(\Lambda(\theta_{1},t))\Lambda_{t}(\theta_{1},t)$. If moreover, at
$\theta_{S}$ there is a split of the agents among different choices, and the
corresponding split fraction is $\lambda_{2}^{(\theta_{S},1,1)}$ (see (11)),
then the entering flow in $B_{0,1}$ through $\theta_{1}$ is $t\mapsto
g(\Lambda(\theta_{1},t))\lambda_{2}^{(\theta_{S},1,1)}(t)\Lambda_{t}(\theta_{1},t)$.
And this is exactly the flow denoted by $g_{01}$ in (10). Similar
considerations (with different function $\Lambda$ and $\lambda$) hold in the
case of agents moving towards $\theta_{2}$ in the branch $B_{1,1}$ and for all
other cases in the other branches.

The above argument is applied in the fixed point procedure in Section 3 for
$\tilde{\psi}_{\varepsilon}$ where, starting from $\rho\in X$ we obtain
another $\rho^{\prime}\in X$, passing through the flow functions. For example
the flow entering in the branch $B_{0,0}$ through $(\theta_{2},0,1)$, that is
the flow $g_{12}$ in (10), using the notation in (12), is (for the
corresponding function $\Lambda$ of the motion from $\theta_{1}$ to
$\theta_{2}$ on the branch $B_{0,1}$) $t\mapsto
g_{0,1}(\Lambda(\theta_{2},t))(\overline{\lambda}_{\varepsilon})_{2}^{(\theta_{1},0,1)}(t)\Lambda_{t}(\theta_{2},t)$,
where $g_{0,1}$ is evaluated as above.


# Anomalous isotope effects of fulleride superconductors

Wei Fan
Key Laboratory of Materials Physics,
Institute of Solid State Physics, Chinese Academy of Sciences,
230031-Hefei, People’s Republic of China

###### Abstract

The numerical calculations of the standard Eliashberg-Nambu strong coupling
theory and a formula of isotope effect derived in this paper provide direct
evidences that an-harmonic vibrations of lattice enhance isotope effect with
anomalous coefficient $\alpha>1/2$. The results in this paper explain very
well the wide distributed $\alpha$ values for the samples with different
ratios of substitutions of 12C by 13C of fulleride superconductor Rb3C60. The
calculations of isotope effects indicate that the intra-molecule radial modes
have more important contributions to superconductivity than the intra-molecule
tangential modes with higher phonon frequencies

PACS: 74.20.Fg, 74.70.Wz, 71.20.Tx, 81.05.Tp

Keywords: Fulleride superconductor, Anomalous isotope effect, An-harmonic
effect, Eliashberg-Nambu Theory

##  1 Introduction

The discoveries of isotope effects in metallic superconductors confirm the
electron-phonon mechanism as the origin of attractive interaction between
electrons. Anomalous isotope effects such as the negative isotope effects with
$\alpha<0$ and the enhanced isotope effects with $\alpha>1/2$ are found in
many new types of superconductors. Understanding the origin of these anomalous
isotope effects is very important to understand the microscopic mechanism of
superconductivity.

Coefficient of isotope effect $0<\alpha<1/2$ is correct for most conventional
superconductors. Negative isotope effect $\alpha<0$, where Tc increases with
mass M, had been found experimentally in conventional superconductors such as
PdH [1, 2] ($\alpha=-0.25$), organic superconductors [3, 4] and
$\alpha$-Uranium ($\alpha=-2.0$) [5]. The negative isotope effect had also
been found in Sr2RuO4 [6] when 16O atoms were substituted by 18O atoms.

Large positive Oxygen isotope effects $\alpha>1/2$ beyond the BCS value $1/2$
had been found in HTSC material La2-xSrxCuO4 ($\alpha$=0.75) at doping level
x=0.12 [7] and in fullerides Rb3C60 ($\alpha>1.0$) with different ratios of
substitutions of 12C with 13C [8, 9, 10]. For La2-xSrxCuO4, near x=0.12 the
large coefficient of isotope effect is probably due to the strong an-harmonic
vibrations companying with structural transitions. Similarly, the transition
from partial substitution to completed substitution of 12C by 13C for Rb3C60
leads to the decrease of $\alpha$ from 1.275 or 1.189 [8, 9] to 0.30 [11]
because the an-harmonic effect is relatively weak at completed substitution.
The Eliashberg theory combined with different models of an-harmonic lattice
vibrations had been used to qualitatively explain these anomalous isotope
effects [12, 13, 14]. Beyond Migdal’s theorem, the coefficient of isotope
effect $\alpha>1/2$ had also been obtained by including non-adiabatic effect
[15].

It’s very convenient to find a formula for isotope effect which is able to
explain anomalous isotope effects such as the enhanced isotope effects with
$\alpha>1/2$ and large negative isotope effect. Most of previous formulas of
isotope effects don’t include the frequency shift due to the isotope
substitution. A formula for isotope effect is derived in this paper using
McMillan’s Tc formula, which explicitly includes the frequency shift of phonon
and is qualitatively consistent with the more accurate numerical results of
Eliashberg-Nambu theory. The formula is similar to a previous formula derived
with similar idea [16]. In this work, we not only compare the results of the
formula with the strong-coupling theory, but also we study in details the
meaning of the an-harmonic coefficient $A^{*}$ and its relation with the
potential $V(r)$ of the vibrations of atoms around their equilibrium positions
in a crystal. We concentrate on the anomalous isotope effects of fulleride
superconductors [17] in the present paper.

An important relation is the key to understand isotope effect and written as

$M\langle\omega^{2}\rangle\lambda=\eta=const$ (1)

where the McMillan-Hopfield parameter $\eta$ characterizes the chemical
environment of atoms in a material [18, 19] and is expressed as
$\eta=N(0)\langle J^{2}\rangle$, where $N(0)$ is the density of state at Fermi
energy, $J$ the matrix element of electron-phonon interaction. Certainly, we
assume $\eta$ is constant for isotope substitution or others small structural
changes. We will see in section (4) that, under the constraint, the an-
harmonic effect is equivalent to the fact that the parameter $\lambda$ is
dependent on $M$ the mass of atom. In a previous work, using electron-phonon
mechanism and the constraint Eq.(1), we have successfully explained the
spatial anti-correlation of the energy gap and the phonon energy of Bi2212
superconductors observed in STM experiments [20]. In this work, using the
constraint combining with McMillan’s Tc formula, we derive a formula of
isotope effect that can explain almost all anomalous isotope effects,
especially for fulleride superconductors.

##  2  Theory

The energy-gap equation of Eliasgberg-Nambu theory in the Matsubara’s
imaginary energy form is standard [18]. When temperature is very close to
transition temperature Tc, the Eliashberg equation can be simplified as

$\displaystyle\sum_{n=0}^{N}(K_{mn}-\rho\delta_{mn})\bar{\Delta}_{n}=0~{}~{}(m\geq
0)$ (2) $\displaystyle K_{mn}$ $\displaystyle=$
$\displaystyle\lambda(m-n)+\lambda(m+n+1)-2\mu^{*}(N)$ $\displaystyle-$
$\displaystyle\delta_{mn}(2m+1+\lambda(0)+2\sum_{l=1}^{m}\lambda(l)),$

where $\bar{\Delta}_{n}=\Delta_{n}/|\omega_{n}|$ is the energy-gap parameter
and
$\lambda(n)=2\int_{0}^{\infty}d\omega\alpha^{2}F(\omega)\omega/(\omega^{2}+(2\pi
nT)^{2})$. The pair-breaking parameter $\rho$ is introduced to form a
eigenvalue problem with $\rho$=0 corresponding the physical energy-gap
equation. The transition temperature Tc is defined as the temperature when the
maximum of eigenvalue of kernel matrix $K_{mn}$ crosses to zero and changes
its sign. We use about N=200 Matsubara energies to solve above equation.

The Eliashberg function is expressed as

$\alpha^{2}F(\omega)=\left\\{\begin{tabular}[]{cc}$\frac{c}{(\omega-\Omega_{P})^{2}+(\omega_{2})^{2}}-\frac{c}{(\omega_{3})^{2}+(\omega_{2})^{2}}$,&$|\omega-\Omega_{P}|<\omega_{3}$\\\
0&others,\end{tabular}\right.$ (3)

where $\Omega_{P}$ is the energy (or frequency) of phonon mode, $\omega_{2}$
the half-width of peak of phonon mode and $\omega_{3}=2\omega_{2}$. We can
write the parameter of electron-phonon interaction
$\lambda=\lambda(0)=2\int_{0}^{\infty}d\omega\alpha^{2}F(\omega)/\omega$. The
moments $\langle\omega^{n}\rangle$ of the distribution function
$(2/\lambda)\alpha^{2}F(\omega)/\omega$ are defined as
$\langle\omega^{n}\rangle=2/\lambda\int^{\infty}_{0}d\omega\alpha^{2}F(\omega)\omega^{n-1}$.
The parameter $\lambda$ characterizes the strength of electron-phonon
interaction by $\lambda\propto N(0)\langle
J^{2}\rangle/M\langle\omega^{2}\rangle$ in terms of Eq.(1). The Coulomb
pseudo-potential is defined as $\mu_{0}=N(0)U$ and its renormalized value as
$\mu^{*}=\mu_{0}/(1+\mu_{0}\ln(E_{C}/\omega_{0}))$, where $U$ is the Coulomb
parameter, $E_{C}$ the characteristic energy for electrons such as the Fermi
energy or band width, and $\omega_{0}$ the characteristic phonon energy such
as the energy cutoff of phonon energy or Debye energy. If
$\omega_{2}\ll\Omega_{P}$, Eq.(1) can be simplified as

$M\Omega_{P}^{2}\lambda=\eta=const,$ (4)

which is obviously correct for Einstein model
$\alpha^{2}F(\omega)=(\lambda/2)\Omega_{P}\delta(\omega-\Omega_{P})$.

##  3  A formula of coefficient of isotope effect

We will derive a formula of isotope effect by including an-harmonic effect.
The shift of phonon energy (or frequency) is explicitly included in the
formula. For single-frequency like mode, the $\alpha$ values in terms of the
formula are qualitatively consistent with the numerical solutions of
Eliashberg equation. Most importantly, the formula clearly shows that the an-
harmonic vibration of lattice leads to the enhanced isotope effect with
$\alpha>1/2$. It’s conveniently to define a parameter $A^{*}$ by
$\Omega_{P}\propto M^{-(1-A^{*})/2}$ to measure the an-harmonic effect. We can
easily see that $A^{*}\neq 0$ represents the an-harmonic effect. The an-
harmonic parameter $A^{*}$ can be expressed as

$A^{*}=1+2M/\Omega_{P}(\delta\Omega_{P}/\delta M).$ (5)

Under harmonic approximation $\Omega_{P}\propto M^{-1/2}$, the parameter
$\lambda$ is a constant ($\delta\lambda/\delta M=0$) because

$A^{*}=1+2M/\Omega_{P}(\delta\Omega_{P}/\delta
M)=-M/\lambda(\delta\lambda/\delta M)=0$

by using the relation $M\Omega_{P}^{2}\lambda=const$.

We start from the McMillan formula of transition temperature Tc of
superconductor

$T_{c}=\frac{\Theta_{D}}{1.45}\exp[-\frac{1.04(1+\lambda)}{\lambda-\mu^{*}(1+0.62\lambda)}].$
(6)

where $\lambda$ is the strength of electron-phonon interaction. The
coefficient $\alpha$ of isotope effect defined by T${}_{c}\propto M^{-\alpha}$
can be obtained from above McMillan formula by the direct mass-dependent from
$\Theta_{D}\propto M^{-1/2}$ and implicit mass-dependent from $\mu^{*}$ by
$\omega_{0}\propto M^{-1/2}$. If an-harmonic effect is included then
$\Theta_{D},\omega_{0}\propto M^{-(1-A^{*})/2}$. In this work we consider
additionally mass-dependent from $\lambda$ by the well known constraint
$M\Omega^{2}_{P}\lambda=\eta$, we obtain a formula of isotope effect which is
expressed as

$\alpha=\frac{1}{2}-\frac{1.04(1+\lambda)(1+0.62\lambda)(\mu^{*})^{2}}{2[\lambda-\mu^{*}(1+0.62\lambda)]^{2}}+A^{*}T(\lambda,\mu^{*})$
(7)

where $A^{*}$ is defined in Eq.(5) and the function $T(\lambda,\mu^{*})$ is
expressed as

$\displaystyle
T(\lambda,\mu^{*})=\frac{1.04\lambda(1+0.38\mu^{*})}{[\lambda-\mu^{*}(1+0.62\lambda)]^{2}}$
(8)

or

$\displaystyle T(\lambda,\mu^{*})$ $\displaystyle=$
$\displaystyle\frac{\lambda(2.08-\lambda)+\lambda(2.7904+1.24\lambda)\mu^{*}}{2[\lambda-\mu^{*}(1+0.62\lambda)]^{2}}$
$\displaystyle+$
$\displaystyle\frac{(0.04+0.42\lambda)(1+0.62\lambda)(\mu^{*})^{2}}{2[\lambda-\mu^{*}(1+0.62\lambda)]^{2}},$

dependent on whether $\Theta_{D},\omega_{0}\propto M^{-1/2}$ for Eq.(8) or
$\Theta_{D},\omega_{0}\propto M^{-(1-A^{*})/2}$ for Eq.(3). We notice that,
for the Eq.(8), the an-harmonic effect enters into the coefficient $\alpha$
only by the M-dependent $\lambda$, however for the Eq.(3) not only by
M-dependent $\lambda$ but also by the $\Theta_{D}$ and $\omega_{0}$. It’s very
important that $T(\lambda,\mu^{*})>0$ for reasonable values
$0<\lambda<2<2.08$, so the sign of the third term of Eq.(7) is determined by
the an-harmonic parameter $A^{*}$. It’s obviously that $A^{*}=0$ under
harmonic approximation $\Omega_{P}\propto M^{-1/2}$ and the third term in the
Eq.(7) is also zero. The derivative $\delta\Omega_{P}/\delta M$ used in the
calculation of an-harmonic parameter $A^{*}$ is approximately obtained from
the experimental energy (or frequency) shift of phonon. The first two terms
give $\alpha<1/2$. The properties of the third term to isotope effect is
determined by the sign of the an-harmonic parameter $A^{*}$. Generally,
$\Omega_{P}$ decreases with increasing mass M, thus
$2(M/\Omega_{P})\delta\Omega_{P}/\delta M<0$. If
$2|(M/\Omega_{P})\delta\Omega_{P}/\delta M|<1$, $A^{*}>0$. So we can get
$\alpha>1/2$ if the third term has larger absolute value than the second term.

##  4  The anomalous isotope effects with $\alpha>1/2$

In many literatures of AnC60, the average mass $m$ of C60 is used in the
definition of coefficient of isotope effect by
$\tilde{\alpha}=-(m/T_{c})\delta T_{c}/\delta m$. If the ratio of substitution
is $p$, the average mass of C60 molecule $m=60[M_{12}(1-p)+M_{13}p]$, so
$\delta m=m-60M_{12}=60p\delta M$ with $\delta M=M_{13}-M_{12}$. The relation
$\tilde{\alpha}=\alpha/p$ connects the parameter $\tilde{\alpha}$ with the
usual definition of $\alpha=-(M/T_{c})(\delta T_{c}/\delta M)$. The usual
coefficient of isotope effect $\alpha=p\tilde{\alpha}=1.275$ is obtained with
$\tilde{\alpha}$=2.125 and $p$=0.60 [8], $\alpha=1.189$ with
$\tilde{\alpha}$=1.45 and $p$=0.82 [9], $\alpha=0.462$ with
$\tilde{\alpha}$=1.4 and $p$=0.33 [10]. For completed substitution with
$p$=1.0, $\alpha=\tilde{\alpha}$=0.30 [11].

The relation $M\langle\omega^{2}\rangle\lambda=\eta=const$ services as a
constraint to determine the parameters $\lambda$ and $\Omega_{P}$ in numerical
calculations. The Coulomb pseudo-potential $\mu^{*}$ has to change in isotope
substitution because $\mu^{*}$ is dependent on the cutoff of phonon energy
$\omega_{0}\propto M^{-1/2}$ or $\propto M^{-(1-A^{*})/2}$ in an-harmonic
approximation. In harmonic approximation, $\delta\mu^{*}=-(\mu^{*})^{2}\delta
M/2M$, the an-harmonic effect is included in calculations only by the
M-dependent $\lambda$. The an-harmonic effect can be realized by shifting
phonon energy $\delta\Omega_{P}$ to make
$A^{*}=1+2M/\Omega_{P}(\delta\Omega_{P}/\delta M)\neq 0$.

The intra-molecule radial mode, which is about $\Omega_{P}$=65 meV or 525 cm-1
from infrared spectrum [10, 8] and the intra-molecule tangential modes around
1400-1 or 174 meV [17] have strong intensity. However we concentrate
attentions on the intra-molecule radial mode with energy $\Omega_{P}$=65 meV
and half-width $\omega_{2}$=8 meV. The Coulomb parameter $\mu_{0}=UN(0)$ and
the corresponding renormalized Coulomb parameter
$\mu^{*}=\mu_{0}/[1+\mu_{0}\ln(E_{C}/\omega_{0})]$ can be estimated from ab-
initio density functional theory (DFT) based on pseudo-potential method using
atomic orbital basis functions[21]. In the DFT calculation, the super-cell
includes one C60 molecule and three Rb atoms. The effect of orientation of
different C60 molecules is ignored. The parameter $U$ is the charge energy
defined as $\delta^{2}E_{tot}/\delta n^{2}=E(n+1)+E(n-1)-2E(n)$, $E(n)$ the
total energy of electric-neutral system of n valence electrons, $E(n+1)$ and
$E(n-1)$ the total energies with one negative and one positive charge
respectively. We choose the possible valence-electron configuration
$4p^{6}5s^{1}$ for Rb atoms and $2s^{2}2p^{2}$for Carbon atoms. The core
electrons are presented by Troulier-Martins pseudo-potentials. The electrons
in semicore state 4p of Rb atoms having already treated as valence electrons
have the single-$\zeta$ basis set and all others valence electrons for all
atoms have the split valence double-$\zeta$ plus polarized basis sets. We use
$\Gamma$ point sampling the first Brillouin zone. The exchange-correlation
potential is GGA Perdew-Burke-Ernzerhof type[22] and the spin-polarization
effects are included in the self-consistent calculations. The results in
table(1) show that the Coulomb parameter $\mu_{0}$ is equal to 4.374 and the
renormalized Coulomb parameters $\mu^{*}$ is 0.127 if the maximum of phonon
energy $\omega_{0}$ is 150meV and $E_{C}$=15 eV when all valence electrons are
included. The common value $\mu^{*}$=0.10 is close 0.127 obtained in this
work. We have known the Coulomb parameter $\mu^{*}$ and the phonon energy
$\Omega_{P}$, the parameter $\lambda$ of electron-phonon interaction can be
defined by the experimental transition temperature. From table(1), we can see
that the parameters $\mu^{*}$ have small influence on isotope effects. Below
we present the calculations of coefficient $\alpha$ of isotope effects with
$\mu^{*}$=0.1.

Table 1:  The Coulomb parameter $\mu^{*}$ is calculated when we choose valence-electron configuration (a) $4p^{6}5s^{1}$ for Rb atom and $2s^{2}2p^{2}$for Carbon atom in DFT calculations. The case (b) is corresponding to the general value $\mu^{*}$=0.10. The coefficients $\alpha$ of isotope effects are calculated using the phonon-energy shifts from 65 meV to 62.3 meV after the isotope substitution (A∗=0.003). | U(eV) | N(0)(1/eV) | UN(0) | $\ln(E_{F}/\omega_{0})$ | $\mu^{*}$ | $\lambda$ | Tc(K) | $\alpha$
---|---|---|---|---|---|---|---|---
(a) | 1.991 | 2.197 | 4.374 | 6.9078 | 0.127 | 0.700 | 29.544 | 0.342
(b) |  |  |  |  | 0.100 | 0.670 | 29.820 | 0.358
Table 2:  The simulation parameters and the results with different an-harmonic parameters $A^{*}$. The coefficients $\alpha$ are calculated based on the numerical solutions of Eliashberg equation and $\alpha^{\prime}$ from the formula Eq.(7) and Eq.(8) derived in this paper. 13C | M (u) |  $\Omega_{P}$ (meV) | $\lambda$ |  $T_{c}$ (K) | $A^{*}$ | $\alpha$ | $\alpha^{\prime}$
---|---|---|---|---|---|---|---
(1) | 13 | 64.3 | 0.632 | 25.5 |  0.742 |  1.929$\pm$0.10 |  2.377
(2) | 13 | 63.4 | 0.650 | 27.2 |  0.409 |  1.207$\pm$0.09 |  1.516
(3) | 13 | 62.3 | 0.674 | 29.0 |  0.003 |  0.358$\pm$0.09 |  0.464
(4) | 13 | 62.2 | 0.676 | 29.2 | -0.034 |  0.282$\pm$0.09 |  0.368
(5) | 13 | 62.0 | 0.680 | 29.5 | -0.108 |  0.132$\pm$0.08 |  0.177
(6) | 13 | 61.8 | 0.685 | 29.8 | -0.182 |  -0.02$\pm$0.08 | -0.014
(7) | 13 | 61.6 | 0.689 | 30.3 | -0.255 | -0.165$\pm$0.08 | -0.205
(8) | 13 | 61.2 | 0.698 | 30.9 | -0.403 | -0.457$\pm$0.08 | -0.588
(9) | 13 | 60.0 | 0.727 | 33.2 | -0.846 | -1.305$\pm$0.08 | -1.736
12C | 12 | 65.0 | 0.670 | 29.82 |  |  |

The parameter of electron-phonon interaction $\lambda$=0.67 is defined by the
experimental transition temperature Tc about 29.5(K) using parameters
$\mu^{*}$=0.1 and $\Omega_{P}$=65 meV. We calculate the $\eta=\eta_{12}$ for
12C. The phonon energies after 13C substitutions are unknown because the
phonon energies are dependent on the ratios of 13C substitutions. We choose
the nine possible energies for phonons in table (2), which are all smaller
than $\Omega_{P}$=65 meV because the energies (or frequency) decrease with
increasing mass of atoms. We also assume that the parameters $\lambda$ of
electron-phonon interaction alter after isotope substitutions. The new
Eliashberg function $\alpha^{2}F(\omega)$ after isotope substitution by 13C is
obtained from the old one before isotope substitution by simply shifting
energy of phonon and scaling it to satisfy the constraint
$M\langle\omega^{2}\rangle\lambda=\eta_{13}=\eta_{12}$. Based on the new
Eliashberg function $\alpha^{2}F(\omega)$, we can calculate the Tc and
$\lambda$ after the isotope substitution. The new Coulomb pseudo-potential
after isotope substitution $\mu^{*}$=0.099583 is obtained from the formula
$\delta\mu^{*}=-(\mu^{*})^{2}\delta M/2M$. The nine possible values of
$\lambda$ corresponding to nine possible energies of phonon are collected in
table (2). The transition temperatures Tc are obtained by solving Eliashberg
equation Eq.(2). The coefficients $\alpha$ are easily calculated in terms of
the transition temperatures Tc after and before isotope substitutions by
$\alpha=-\ln(^{13}T_{c}/^{12}T_{c})/\ln(M_{13}/M_{12})$. The an-harmonic
parameters $A^{*}$ are calculated in terms of the shifts of phonon energies.

For the case (2) in table (2), we get $\alpha=1.207$ which is very close to
the already known maximum value 1.275 in experiments with uncompleted
substitution ($\tilde{\alpha}$=2-2.25, $p$=0.60) [8]. For sample with
completed substitution with 13C, the value of $\alpha$ decreases to $\sim$0.30
[11]. To explain the interesting results, we assume the increase of ratio of
substitution makes the distribution of 13C more homogenous and the an-harmonic
effect become weak with small $A^{*}$. If the phonon energy shifts to 62.3 meV
with smaller an-harmonic parameter $A^{*}$=0.003, the value of $\alpha$=0.358
is close to the experimental values from 0.30 to 0.37 for the samples with
completed substitution.

Figure 1:  The comparison of coefficients of isotope effect having obtained
from numerical calculations of Eliashberg equation and from Eq.(7) for radial
intra-molecule mode $\Omega_{P}$=65 meV (a) and the tangential intra-molecule
mode $\Omega_{P}$=150 meV (b)

For negative an-harmonic parameter$A^{*}$=-0.182, we find $\alpha$=0.02 close
to zero. After having made $A^{*}$ more negative further, at $A^{*}$=-0.403,
we get the negative coefficient $\alpha$=-0.457. However, the negative isotope
effect wasn’t found in fulleride superconductors. From above calculations we
can see that the sign of $A^{*}$ determines the sign of $\alpha$ at relative
larger absolute values of $A^{*}$. If the absolute values of an-harmonic
parameters are not too large $0.07>A^{*}>-0.21$ or an-harmonic effect is weak,
the coefficient $\alpha$ are from 0 to 0.5 within the range of general values.
The above calculations explain the wide distributed values of coefficients for
Rb3C60 in experiments. This is because the an-harmonic parameters are
dependent on the ratios of substitutions.

From table (2), we can see that the values obtained from numerical
calculations are very close to the values from Eq.(7) with Eq.(8). The values
of $\alpha$ and $\alpha$’ have the same sign on the same rows. It’s very
important that if $\alpha>1/2$, $A^{*}>0$ must be satisfied. Thus, the
enhanced coefficients $\alpha$ larger than $1/2$ for Rb3C60 in uncompleted
substitution samples mean that the strong an-harmonic effects. If the an-
harmonic effect becomes weak as approaching the completed substitution, the
value of $\alpha$ will decrease and reach to smaller value about 0.358 at
$A^{*}=0.003$.

Fig. 1(a) shows the results of the numerical solutions of Eliashberg equation
and the formula Eq.(7). We can see that Eq.(7) with Eq.(8) is more close to
the numerical solution than using the Eq.(3) although for Eq.(3) more
completed an-harmonic effects are included. We have preformed the same
calculations for intra-molecule tangential mode with energy $\Omega_{P}$=150
meV, $\mu^{*}$=0.285, $\lambda$=0.67 and Tc=29.5(K). From Fig.1(b), there is
larger slope of $\alpha$-$\delta\Omega_{P}/\Omega_{P}$ curve compared with the
intra-molecule radial mode with energy $\Omega_{P}$=65 meV. The coefficient of
isotope effect of the intra-molecular radial mode $\alpha$=0.358 at
$A^{*}\simeq$0.0 is more close to experiments at completed substitutions
$\alpha\sim 0.30-0.37$. So it is more correlated with superconductivity of
fulleride than the intra-molecule tangential mode.

##  5 An-harmonic parameter $A^{*}$ and models of lattice vibrations

Figure 2:  An-harmonic coefficients $A^{*}$ for the models of interaction with
form $V(r)=hr^{d}$. The inserted figures (1),(2) and (3) illustrate the shapes
of interacting potentials and the quantized energy levels for $h$=0.01 and
$d$=0.5, 2.0 and 4.0 respectively. The maximum of quantum number of angular
momentum is up to $l_{max}$=8

In this paper, we introduce an important parameter, the an-harmonic parameter
$A^{*}$. To make parameter $A^{*}$ more realistic, we study the models of
lattice vibrations in a crystal. The atoms in crystal move around their
equilibrium positions. The change of potential energy of an atom away from its
equilibrium position is written as $V(r)=hr^{d}$ with radially local
displacement $r$. The harmonic vibration is corresponding $d=2$. We had
numerically solved the Schrödinger equation
$[-\hbar^{2}\nabla^{2}/2M+V(r)]\psi(\vec{r})=E\psi(\vec{r})$ under isotropic
approximation, especially, the one-dimensional radial part converts into
difference equation. We calculated the relation $E_{0}(M)$ between energy of
ground state and mass of atom for a series of discrete values of $M$. It’s
well known that for harmonic approximation $E_{0}(M)\propto 1/M^{0.5}$. The
effective an-harmonic parameters $A^{*}$ for interaction with the form
$V(r)=hr^{d}$ are obtained by fitting the discrete $E_{0}(M)$ functions with
function $\tilde{E}_{0}(M)\propto 1/M^{(1-A^{*})/2}$. From the Fig.2, if $d<2$
we can obtain $A^{*}>0$. So the enhanced isotope effects with $\alpha>0.5$ are
hopefully found based on the formula Eq.(7) obtained in this work. If $d>2$ so
$A^{*}<0$, we can get the normal isotope effects with $\alpha<0.5$ and the
negative isotope effects with $\alpha<0$. We can see that the simple model
$V(r)=hr^{d}$ is suitable for the isotope effects of fulleride Rb3C60 for
different ratios of substitutions because the an-harmonic parameters are
within the range to obtain experimental values of $\alpha$.

##  6  Discussion and Summary

The Figure 1(b) shows clearly that the formula Eq.(7) isn’t good approximation
to more accurate numerical solution for high-energy mode. The reason is
probably that McMillan’s formula isn’t correct in some regions of parameter
space. The McMillan’s formula is good approximation when the parameter
$\lambda$ of electron-phonon interaction is not too large ($\lambda<1$) and
the phonon energy $\Omega_{P}$ is not too high. The inter-molecule phonon
modes such as the vibration between $C_{60}$ molecules and the between alkali-
metal atoms and $C_{60}$ molecules are ignored in this work because their
energies generally smaller than 12 meV. To obtain $T_{c}$=30K, the parameters
$\lambda$ of electron-phonon interaction are at least 3.0 so the instability
of lattice will destroy superconductivity. However, these inter-molecule modes
still have significant influence on the properties of fulleride
superconductors such as the differences of isotope effects for different
substitution configurations $Rb_{3}[(^{13}C_{60})_{x}(^{12}C_{60})_{1-x}]$ and
$Rb_{3}(^{13}C_{x}^{12}C_{1-x})_{60}$ [11].

In summary, the coefficients $\alpha$ of isotope effects obtained in this
paper are very close to the values in experiments or within the range of
experimental values. The reductions of $\alpha$ with increasing the
substitution ratios of 13C are due to the reductions of an-harmonic effects of
lattice vibrations when the ratios of substitution tend to 100%. The enhanced
coefficients of isotope effects with $\alpha>1/2$ generally happen at the
intermediate stage of transition from one phase to another. Finally, the
formula Eq.(7) and the numerical methods used in this paper are also suitable
to study isotope effects of other superconductors if the electron-phonon
interaction is the pairing mechanism for electrons.

The author benefits from the discussions with members within CMS Group at
ISSP-CAS, especially from Dr. Xu Yong, Huang Ling-Feng, Li Long-Long and Li
Yan-Ling. This work is supported by Director Grants of Hefei Institutes of
Physical Sciences, Knowledge Innovation Program of Chinese Academy of Sciences
and National Science Foundation of China.

## References

  * [1] T. Skoskiewicz, Phys. Status Solidi A11, K123 (1992).
  * [2] R. J. Miller and C. B. Satterthwaite, Phys. Rev. Lett. 34, 144 (1975).
  * [3] J. A. Schlueter, J. M. Williams, A. M. Kini, U. Geiser, J. D. Dudek, M. E. Kelly, J. P. Flynn, D. Naumann, and T. Roy, Physica C265, 163 (1996).
  * [4] A. M. Kini, K. D. Carlson, H. H. Wang, J. A. Schlueter, J. D. Dudek, S. A. Sirchio, U. Geiser, K. R. Lykke and J. M. Williams, Physica C264, 81 (1996).
  * [5] R. D. Fowler, J. D. G. Lindsay, R. W. White, H. H. Hill, and B. T. Matthias, Phys. Rev. Lett. 19, 892 (1967).
  * [6] Z. Q. Mao, Y. Maeno, Y. Mori, S. Sakita, S. Nimori, and M. Udagawa, Phys. Rev. B63, 144514 (2001).
  * [7] M. K. Crawford, M. N. Kunchur, W. E. Farneth, E. M. McCarron III and S. J. Poon, Phys. Rev. B41, 282 (1990).
  * [8] A. A. Zakhidov, K. Imaeda, D. M. Petty, K. Yakushi, H. Inokuchi, K. Kikuchi, I. Ikemoto, S. Suzuki, and Y. Achiba, Phys. Lett. A164 , 355 (1992).
  * [9] P. Auban-Senzier, G. Quirion, D. Jerome, P. Bennier P, S. Della-Negra, C. Fabre and A. Rassat, Synth. Met. 56, 3027 (1993).
  * [10] T. W. Ebbesen, J. S. Tsai, K. Tanigaki, J. Tabuchi J, Y. Shimakawa, Y. Kubo, I. Hirosawa and J. Mizuki, Nature 355, 620 (1992).
  * [11] C. C. Chen, and C. M. Lieber, Science 259, 655 (1993).
  * [12] T. Galbaatar, S. L. Drechsler, N. M. Plakida, and G. M. Vujicic, Physica C176, 496 (1991).
  * [13] T. Galbaatar, R. -J. Rakauska, and J. Sulskus, International Journal of Modern Physics B3, 897 (1989).
  * [14] Y. Inada and K. Nasu, J. Phys. Soc. Jpn. 61, 4511 (1992).
  * [15] C. Grimaldi, L. Pietronero, and S. Strässler, Phys. Rev. B52, 10530 (1995).
  * [16] V. H. Crespi, M. L. Cohen, D. R. Penn, Phys. Rev. B43, 12921 (1991).
  * [17] O. Gunnarsson, Reviews of Modern Physics 69, 575 (1997).
  * [18] P. B. Allen and R. C. Dynes, Phys. Rev. B12, 905 (1975).
  * [19] W. L. McMillan, Phys. Rev. 167, 331 (1968).
  * [20] W. Fan, Chinese Physics Letter, 25, 2217 (2008).
  * [21] J. M. Soler, E. Artacho, J. D. Gale, A. García, J. Junquera, P. Ordejón and D. Sánchez-Portal, J. Phys. Condens. matter 12, 2745 (2002).
  * [22] J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).