<figure><img src="image_1.png" /></figure>

HAL
open science

# Defective Laplacians and paradoxical phenomena in crowd motion modeling

Bertrand Antti Maury

<figure><img src="image_2.png" /></figure>

► To cite this version:

Bertrand Antti Maury. Defective Laplacians and paradoxical phenomena in crowd motion modeling. Comptes Rendus. Mécanique, In press, Comptes Rendus. Mécanique. Académie des Sciences, 10.5802/crmeca.205. hal-04138286

HAL Id: hal-04138286

https://hal.science/hal-04138286

Submitted on 22 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Defective Laplacians and paradoxical phenomena in crowd motion modeling

Bertrand Maury<sup>a</sup>

<sup>a</sup> LMO, Université Paris-Saclay, F-91405 Orsay cedex & DMA, Ecole Normale Supérieure, PSL University, Paris

E-mail: bertrand.maury@universite-paris-saclay.fr (B. Maury)

**Abstract.** In both continuous and discrete settings, Laplace operators appear quite commonly in the modeling of natural phenomena, in several contexts: diffusion, heat propagation, porous media, fluid flows through pipes, electricity. ... In these contexts, the discrete Laplace operator enjoys all the properties of its continuous counterpart, in particular: self-adjointness, variational formulation, stochastic interpretation, mean value property, maximum principle, ... In a first part, we give a detailed description of the correspondence between these mathematical properties and modeling considerations, in contexts where the continuous and the discrete settings perfectly match. In a second part, we describe a pathological situation, in the context of granular crowd motion models. Accounting for the non-overlapping constraint between hard discs leads to a particular operator acting on a field of Lagrange multipliers, defined on the dual graph of the contact network. This operator is defective in a certain sense: although it is the microscopic counterpart of the macroscopic Laplace operator, this discrete operator indeed lacks some properties, in particular the maximum principle. We investigate here how this very *defectivity* may explain some *paradoxical* phenomena that are observed in crowd motions and granular materials, phenomena that are not reproduced by macroscopic models.

**Résumé.** Aux niveaux continu et discret, l'opérateur de Laplace intervient de façon très courante dans la modélisation de phénomènes naturels, dans de nombreux contextes: diffusion, propagation de la chaleur, milieux poreux, écoulement de fluides dans des conduits, électricité. ... Dans ces contextes, le laplacien discret possède toutes les propriétés de son pendant continu: caractère auto-adjoint, structure variationnelle, interprétation stochastique, propriété de la valeur moyenne, principe du maximum, .... Dans une première partie, nous proposons une description détaillée des liens entre ces propriétés et les aspects de modélisations, dans des contextes où les notion continues et discrètes se correspondent parfaitement. Dans une seconde partie, nous décrivons une situation pathologique, dans le contexte de la modélisation de foules d'un point de vue granulaire. La prise en compte de la contrainte de non recouvrement entre grains rigides conduit à un opérateur particulier qui agit sur les champs de multiplicateurs de Lagrange, définis sur le graphe dual du réseau de contacts. Cet opérateur est déficient dans un certain sens : bien qu'il apparaisse comme le pendant discret du laplacien continu, il ne vérifie pas certaines des propriétés usuelles du laplacien, en particulier le principe du maximum. Nous explorons comment cette déficience permet d'expliquer certains effets paradoxaux observés en mouvements de foules, que les modèles continus ne reproduisent pas.

**Keywords.** Discrete Laplace operator, maximum principle, crowd motion, faster-is-slower effect.

**Mots-clés.** Laplacien discret, principe du maximum, mouvements de foules, effet faster-is-slower.
# 1. Introduction, continuous and discrete Laplacians

The Laplace operator occupies a central position in the modeling of physical phenomena. In many situations, it combines a phenomenological phenomenon together with a conservation principle. In the context of diffusion of some substance represented by its density $p = p(x, t)$, Fick's Law states that the flux vector $J$ is proportional to the gradient of $p$:

$$J = -D\nabla p,$$

where $D > 0$ is the diffusion coefficient. The flux vector $J$ is such that, for any oriented surface $\Sigma$, with normal vector $n$, $\int_{\Sigma} J \cdot n$ is the flux of substance through $\Sigma$, counted positively in the $n$-direction. Now considering a domain $\omega$, assuming that there is a source term $f$, mass balance writes

$$\frac{d}{dt} \int_{\omega} p = \int_{\omega} \partial_t p = \int_{\omega} f - \int_{\omega} J \cdot n = \int_{\omega} f - \int_{\omega} \nabla \cdot J$$

which leads to

$$\int_{\omega} (\partial_t p + \nabla \cdot J - f) = 0.$$

Since $\omega$ is arbitrary, and $J = -D\nabla p$, it yields the heat equation

$$\partial_t p - D\Delta p = f,$$

the stationary version of which is the Poisson equation

$$-D\Delta p = f,$$

set in a domain $\Omega$ occupied by the substance. Let us mention another context in which $-\Delta$ appears as a composition of $\nabla$ (gradient) and $-\nabla \cdot$ (minus the divergence operator). We consider the flow of an incompressible viscous fluid in a porous medium. Darcy's law writes

$$u = -k\nabla p$$

where $p$ is the pressure, $k$ the permeability of the medium, and $u$ is the Darcy velocity (which should be thought of as a volume flux, expressed in cubic meter per square meter per second, which is indeed homogeneous to a velocity). Volume conservation writes

$$-\nabla \cdot u = -\nabla \cdot k\nabla p = 0,$$

which is again the Laplace equation in a possibly non homogeneous medium.

To obtain a discrete Laplacian, one may e.g. consider a set of cells filled with some medium (see Fig. 1). A substance freely diffuses in each cell, so that the concentration within a cell can be considered as uniform. Cells are separated by a membrane, and the flux of substance between two cells is assumed to be proportional to the difference between concentrations. We define the abstract non-oriented graph $(V, E)$ as follows: the set of vertices $V$ is the set of cells (symbolized by black dots in the figure), and the set of edges $E$ contains $(x, y)$ if and only if cells $x$ and $y$ share a piece of membrane. Note that $E$ is symmetric by construction: the graph is said to be *non oriented*. We shall write $x \sim y$ whenever $(x, y) \in E$. If we denote by $p_x$ the concentration (or partial pressure) on the cell $x$, and $c_{xy}$ the permeability of the membrane ($c_{xy}$ plays the role of a *conductance* in electrical / fluid networks), the flux from $y$ to $x$ writes

$$J_{y \to x} = c_{xy}(p_y - p_x).$$

Note that the flux is symbolically associated to the edge connecting $x$ and $y$, but it corresponds here to a transfer that is spread over the interface which separates $x$ and $y$ cells (grey arrows in the figure).
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Diffusion process between cells.</figcaption></figure>

Denoting by $M_x$ the volume of cell $x$, the product $M_x p_x$ is the total quantity of substance at $x$, and mass balance at each cell $x$ writes

$$M_x \frac{dp_x}{dt} = - \sum_{y \sim x} c_{xy}(p_x - p_y).$$

It can be written globally

$$M \frac{dp}{dt} + Lp = 0,$$

where $L$ is a discrete Laplacian, i.e. the discrete counterpart of $-\Delta$, defined as

$$L : p \in \mathbb{R}^V \to Lp = \left( \sum_{y \sim x} c_{xy}(p_x - p_y) \right)_{x \in V},$$

and $M = \text{diag}(M_x)$.

In another context, natively static, one may also consider a network of interconnected pipes, in which a viscous fluid is flowing. If we denote by $V$ the set of bifurcation points, by $E \subset V \times V$ the (symmetric) set of pipes between vertices, the flow through a pipe $(x, y)$ follows the so-called Poiseuille's law

$$u_{xy} = c_{xy}(p_x - p_y),$$

where $c_{xy}$ is the conductance of the edge $(x, y)$, which depends on geometrical properties of the actual pipe that it represents. Now consider that the flow in conservative at some vertices, the set of which is denoted by $\mathring{V}$, and possibly non conservative on the set of remaining points $\Gamma = V \setminus \mathring{V}$ (see Fig. 2). Conservation on $\mathring{V}$ writes

$$\sum_{y \sim x} u_{xy} = \sum_{y \sim x} c_{xy}(p_x - p_y) = 0 \quad \forall x \in \mathring{V}.$$

As a consequence, if the value of the pressure is prescribed at boundary points (field $P \in \mathbb{R}^\Gamma$), the pressure field is the solution to a Dirichlet problem

$$\left| \begin{array}{ll} Lp = 0 & \text{in } \mathring{V} \\ p = P & \text{on } \Gamma. \end{array} \right.$$
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Network $(V, E, c, \Gamma)$</figcaption></figure>

This setting has been intensively used to model the airflow in the respiratory tree (see e.g. [1] or [2]). It can be transposed to resistive electrical networks, where the potential plays the role of the pressure, and electric intensity the role of fluid fluxes.

**Remark 1.** Let us make it clear that the discrete equations written above do not aim here at approximating the solution to continuous Poisson or Laplace equations, they model intrinsically discrete phenomena, and so will the discrete Laplace operator introduced in Section 3. Let us nevertheless mention the deep link between the diffusion between cells which we considered (Fig. 1) and the elaboration of numerical methods of the Finite Volume type, for which cells are virtual, since they do not correspond in general to a specific physical zone of the underlying domain. We refer to [3] for a very detailed and documented account of discretization strategies of the Finite Volume type. Requiring that the discrete operator verify the maximum is essential to establish $\ell^\infty$ stability of the numerical scheme. Note that, in the context of approximating PDE's, convergence also relies on consistency with respect to the continuous operator, which calls for strong requirements on the mesh that is used. In the present setting, as detailed in the next section, the maximum principle is automatically verified since the models are compliant with the principles of thermodynamics, whereas consistency with respect to the continuous Laplacian is not an issue.

## 2. Mathematical properties and their physical interpretations

We recall in this section the main properties of the standard Laplace operator on euclidean domains, and formulate their discrete counterparts for resistive networks. Some of the results at the continuous level require the domain to be *smooth* in some sense. We shall simply call smooth a domain with Lipschitz boundary.

### 2.1. Self-adjointness, positivity, variational setting

Consider the case of Dirichlet boundary conditions, i.e. $p$ is set to 0 on $\partial\Omega$. From a functional analysis standpoint, $-\Delta$ is an unbounded operator in $L^2(\Omega)$, which is self-adjoint, thanks to Green's formula:

$$ \langle -\Delta p | q \rangle = - \int_{\Omega} q \Delta p = \int_{\Omega} \nabla p \cdot \nabla q = \langle -\Delta q | p \rangle. $$
Its self adjoint character comes from the fact that “$-\nabla$” (opposite of the divergence) and $\nabla$ (gradient) are mutually adjoint for the $L^2$- duality. Note that, although it is intensively used in the theoretical study of linear elliptic equations, this fact is in some way *fortuitous*, and makes a connections between considerations of very different natures. On the one hand, the divergence operator expresses a universal conservation principle, sometimes called Lavoisier' Principle. It would not make any sense to imagine an alternative way, e.g. nonlinear, to express this conservation principle. On the other hand Fick's law is phenomenological, it was inferred from experimental measurements, and in some contexts it takes a nonlinear form (see e.g. [4]).

This makes it possible to formulate the Dirichlet problem as a minimization problem, the Poisson equation corresponding to optimality conditions. It can be expressed in the form of a very classical existence and uniqueness property, which we shall write extensively anyway to emphasize the similarity with the discrete context (see Prop. 5 thereafter).

**Proposition 2.** Let $\Omega$ be a bounded smooth domain, and let $P \in H^{1/2}(\Gamma)$ be given. Then the Dirichlet problem

$$ \left| \begin{array}{ll} -\Delta p = 0 & \text{in } \Omega \\ p = P & \text{on } \Gamma. \end{array} \right. \qquad (1) $$

admits a unique weak solution in $H^1(\Omega)$, in the sense that there exists a unique $p \in H_p^1(\Omega)$ such that

$$ \int_{\Omega} \nabla p \cdot \nabla q = 0 \quad \forall q \in H_0^1(\Omega), $$

with

$$ H_p^1(\Omega) = \{q \in H^1(\Omega), q|_{\Gamma} = P\}, $$

where $q|_{\Gamma}$ denotes the trace of $q$ on $\Gamma$. Furthermore, $p$ is the unique minimizer of

$$ v \mapsto J(v) = \frac{1}{2} \int_{\Omega} |\nabla q|^2 $$

in $H_p^1$.

**Proof.** See e.g. [5]. □

<figure><img src="image_6.png" /></figure>

In the discrete setting, a similar framework can be set. Let us consider a resistive network $(V, E, c, \Gamma)$, where $c \in [0, +\infty]$ is the collection of conductances, and $\Gamma$ is the subset of vertices which may exchange fluid (or intensity in the electrical context) with the outside world. We shall consider pressure fields $p \in \mathbb{R}^V$ and flux fields $u \in \mathbb{R}^E$, with the convention that fluxes are skew symmetric, i.e. $u_{xy} = -u_{yx}$. As detailed in [1], we can define a discrete divergence operator (it actually stands for the *opposite* of the divergence):

$$ \begin{gathered} \partial : u \in \mathbb{R}^E \longrightarrow \partial u \in \mathbb{R}^V \\ (\partial u)_x = \sum_{y \sim x} u_{yx}. \end{gathered} $$

The conservation at interior vertices simply writes $(\partial u)_x = 0$ for any $x \in \tilde{V}$.

In the same spirit, we define a discrete gradient (or “drop”) operator as

$$ \begin{gathered} \partial^* : p \in \mathbb{R}^V \longrightarrow \partial^* p \in \mathbb{R}^E \\ (\partial^* p)_{xy} = p_y - p_x. \end{gathered} $$

Poiseuille / Ohm's law then writes $u = -c\partial^* p$.

The two operators $\partial : \mathbb{R}^E \to \mathbb{R}^V$ and $\partial^* : \mathbb{R}^V \to \mathbb{R}^E$ defined above can be straightforwardly checked to be mutually adjoint, as stated by the next proposition.
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Abstract setting: chains and cochains.</figcaption></figure>

**Proposition 3.** We have that

$$ \langle \partial u | p \rangle_V = \langle u | \partial^* p \rangle_E \quad \forall u \in \mathbb{R}^E, p \in \mathbb{R}^V, $$

where $\langle \cdot | \cdot \rangle_V$ and $\langle \cdot | \cdot \rangle_E$ are the canonic scalar products on $\mathbb{R}^V$ and $\mathbb{R}^E$, respectively.

**Proof.** For $u \in \mathbb{R}^E, p \in \mathbb{R}^V$, it holds that

$$ \langle \partial u | p \rangle_V = \sum_{x \in V} p_x(\partial u)_x = \sum_{x \in V} p_x \sum_{y \sim x} u_{yx} = \sum_{(x,y) \in E} u_{xy} \underbrace{(p_y - p_x)}_{\partial^* p(x,y)} = \langle u | \partial^* p \rangle_E, $$

which ends the proof.

<figure><img src="image_8.png" /></figure>

**Remark 4.** In the language of algebraic topology, the operators $\partial$ and $\partial^*$ can be considered as boundary and coboundary operators, respectively, whereas fluxes are 1-chains, mass conservation defects are 0-chains, pressure fields are 0-cochains, and pressure drops are 1-cochains (see Fig. 3). The duality between 1-chains and 1-cochains (on the left-hand side of the figure) corresponds to the power dissipated within edges, whereas the duality between 0-chains and 0-cochains corresponds to the power of external pressures (i.e. pressure at non-conservative vertices where some fluid is exchanged with the outside world). In the spirit of Remark 1, it is fruitful to establish a link between this abstract setting and discretization strategies for Partial Differential Equations. We refer to [6] for a detailed description of how principles borrowed from exterior calculus can be integrated in the elaboration of Finite Element discretization strategies which are respectful of the underlying equations, from this standpoint.

The problem takes the form of a discrete Darcy problem, which combines a phenomenological law (Poiseuille, Fick, or Ohm), together with a conservation principle:

$$ \left| \begin{array}{ll} u + c\partial^* p = 0 & \text{over } E \\ \partial u = 0 & \text{in } V. \end{array} \right. \qquad (2) $$

Like in the continuous setting, the discrete Laplace operator appears through the elimination of the velocity in the Darcy system. It can be written $L = \partial c \partial^*$ where, for any $p \in \mathbb{R}^V$, $c \partial^* p$ simply stands for the term-by-term product of $c$ and $\partial^* p$, i.e. $c \partial^* p \in \mathbb{R}^E$, with $(c \partial^* p)_e = c_e (\partial^* p)_e$. Its self-adjointness is a straight consequence of Proposition 3:

$$ \langle Lp | q \rangle = \langle \partial c \partial^* p | q \rangle = \langle c \partial^* p | \partial^* q \rangle = \langle \partial^* p | c \partial^* q \rangle = \langle p | \partial c \partial^* q \rangle = \langle p | Lq \rangle. $$

Like in the continuous setting again, this property can be seen as fortuitous in terms of modeling, since it relies on a strong relation between a universal principle (mass conservation) and a phenomenological law (Poiseuille's in a fluid context, or Ohm's law for an electrical network).

Proposition 2 has a discrete counterpart: the Dirichlet problem is well-posed under general assumptions, and the problem can be formulated in a variational way.

**Proposition 5.** Let $(V, E, c, \Gamma)$ be a connected resistive network, with boundary $\Gamma \neq \emptyset$, and $P \in \mathbb{R}^\Gamma$ a collection of boundary value. Then the Dirichlet problem

$$ \left| \begin{array}{l} (Lp)_x = \sum_{y \sim x} c_{xy} (p_x - p_y) = 0 \quad \forall x \in \mathring{V} \\ p_x = P_x \quad \forall x \in \Gamma. \end{array} \right. \qquad (3) $$
admits a unique solution in $\mathbb{R}^V$. It is the unique field $p \in H_P$ such that

$$ \sum_{(x,y) \in E} c_{xy}(p_x - p_y)(q_x - q_y) = 0 \quad \forall q \in H_0, $$

with

$$ H_P = \{q \in \mathbb{R}^V, q|_{\Gamma} = P\}. $$

Furthermore, $p$ is the unique minimizer of

$$ q \mapsto J(q) = \frac{1}{2} \sum_{(x,y) \in E} c_{xy}(q_x - q_y)^2 $$

in $H_P^1$.

**Proof.** Any field $q$ of $\mathbb{R}^V$ can be written $(\hat{q}, q_Γ)$ so that minimizing $J$ over $H_P$ amounts to minimizing

$$ \hat{q} \mapsto \Phi(\hat{q}) = J(\hat{q}, P_Γ) $$

over $\mathbb{R}^{\hat{V}}$. Let us show that $\Phi$ is coercive. Let $q = (\hat{q}, P_Γ)$ be given in $H_P$. Since the graph is connected and $Γ \neq \emptyset$, any $x \in \hat{V}$ can be connected to $Γ$ by a path

$$ x_0 \in Γ \sim x_1 \sim \dots \sim x_n = x. $$

It holds that

$$ \begin{aligned} |q_x| &= |q_{x_0} + \sum_{j=1}^{n} (q_{x_{j-1}} - q_{x_j})| \le |q_{x_0}| + \sum_{j=1}^{n} |q_{x_{j-1}} - q_{x_j}| \\ &\le |P_{x_0}| + \underbrace{\left( \sum_{j=1}^{p} c_{x_{j-1}x_j} |q_{x_{j-1}} - q_{x_j}|^2 \right)^{1/2}}_{\le J(\hat{p})^{1/2}} \left( \sum_{j=1}^{p} \frac{1}{c_{x_{j-1}x_j}} \right)^{1/2}. \end{aligned} $$

Boundedness of $\Phi(\hat{q})$ therefore implies boundedness of $q$ (for any norm in this finite dimensional context). The functional $\Phi$ is therefore coercive: $\Phi(\hat{q}) \to +\infty$ when $|q| \to +\infty$. Since it is continuous, it admits a minimizer $\hat{p}$, at which $\nabla\Phi$ vanishes, so that $p = (\hat{p}, P_Γ)$ solves (3), which ensures the existence of a solution. Let us now prove that $\Phi$ is strictly convex. Considering $\hat{p}$ and $\hat{p}'$ two different fields, we consider $p = (\hat{p}, P_Γ)$ and $p' = (\hat{p}', P_Γ)$. Since $p \neq p'$, while they identify at least at one vertex, and since the graph is connected, there exist at least an edge $(x, y)$ on which $p_x - p_y \neq p'_x - p'_y$. The associated quadratic term of $\Phi$ therefore ensures strict convexity. The minimizer of $\Phi$ is therefore unique. By convexity of $\Phi$, any field $\hat{p}$ that vanishes the gradient (i.e. such that $(\hat{p}, P_Γ)$ is a solution to (3)) is a minimizer, which ensures that the solution to (3) is unique. $\square$

Let us insist on the fact that, in the discrete setting, the notion of boundary is somewhat arbitrary, whereas for a euclidean domain its designs the topological boundary $\bar{\Omega} \setminus \Omega$. The latter expression makes no sense for a network, and $Γ$ is simply the set of non conservative vertices. We may even consider the caricatural case $Γ = V$, $\hat{V} = \emptyset$, which leads to a trivial Dirichlet problem $p \equiv P$.

In both the continuous and the discrete settings, the positivity of the Laplace operator that results from the Darcy system expresses dissipation of energy. For the continuous Darcy problem, $c|\nabla p|^2$ is the dissipated power per unit volume (or area in 2d), due to viscous stress in the fluid. At the discrete level, $(c|\partial^* p|^2)_e$ is the dissipated power within edge $e$, through dissipative stress in the fluid context, or through Joule's effect in the electrical context ($p$ then represents an electrical potential).
## 2.2. Mean value property and maximum principle

The Laplacian of a function at some point measures the discrepancy between the value at this point and the mean value around it, as stated by the following proposition.

**Proposition 6.** Let $\Omega \in \mathbb{R}^d$ be a domain (nonempty open set). A scalar field $p$ is an harmonic function in $\Omega$, i.e. $-\Delta p \equiv 0$ in $\Omega$, if and only if for any $x \in \Omega$, any $r > 0$ such that $\bar{B}(x, r) \subset \Omega$,

$$p(x) = \frac{1}{|B(x, r)|} \int_{B(x, r)} p(y) dy = \frac{1}{|S(x, r)|} \int_{S(x, r)} p(y) dy,$$

where $|B(x, r)|$ is the $d$-volume of the ball $B(x, r)$, and $|S(x, r)|$ the $d-1$ dimensional measure of the sphere $S(x, r) = \partial B(x, r) \in \mathbb{R}^d$.

**Proof.** See [7], Chapter I.

<figure><img src="image_9.png" /></figure>

Beyond this characterization of harmonic functions, the Laplacian operator estimates the local gap between the value of a function at some point $x$ and the mean value in the neighborhood of $x$. This general remark can be formulated as follows:

**Proposition 7.** Let $p$ be a $C^4$ scalar field in an open set containing $x \in \mathbb{R}^d$. Then

$$-\Delta p(x) = \frac{1}{\epsilon^2} \left( p(x) - \frac{1}{|S(x, \epsilon)|} \int_{S(x, \epsilon)} p(y) dy \right) + \mathcal{O}(\epsilon^2).$$

**Proof.** For $\epsilon > 0$ sufficiently small, for any $\sigma \in S(0, 1)$, it holds that

$$f(x + \epsilon\sigma) - f(x) = \sum_{|\alpha|=1}^{3} \frac{\partial^\alpha f}{\partial x^\alpha}(x) \epsilon^{|\alpha|} \sigma^\alpha + \mathcal{O}(\epsilon^4),$$

with $|\alpha| = \alpha_1 + \cdots + \alpha_d$ for $\alpha = (\alpha_1, \dots, \alpha_d) \in \mathbb{N}^d$, and $\sigma^\alpha = \sigma_1^{\alpha_1} \times \dots \sigma_d^{\alpha_d}$. We integrate over $S(0, 1)$, and we use the fact that the terms associated to $|\alpha| = 1$, $|\alpha| = 3$, and $|\alpha| = 2$ except for the diagonal terms (trace of the Hessian matrix), all vanish by symmetry of $S(0, 1)$. The proposed expansion of $-\Delta p(x)$ is obtained by dividing the Taylor expansion by the measure of $S(x, \epsilon)$.

<figure><img src="image_10.png" /></figure>

**Proposition 8.** Let $\Omega$ be a smooth domain, and $p$ harmonic in $\Omega$. Then the maximum of $p$ is attained on the boundary of $\Omega$.

**Proof.** This is a straight consequence of Prop. 6.

<figure><img src="image_11.png" /></figure>

Another type of maximum principle, namely Hopf maximum principle, will be important to interpret the behaviour of macroscopic crowd motion models

**Proposition 9.** Let $\Omega$ be a smooth domain, and $p$ a smooth scalar field such that $-\Delta p \ge 0$. At any point of $\partial\Omega$ which minimizes $p$, it holds that $\partial p/\partial n < 0$.

**Corollary 10.** Let $\Omega$ be a smooth domain, and $p$ such that $-\Delta p = f$, with $f$ a smooth nonnegative function that is non identically zero. Then $p$ is positive in $\Omega$, and $\partial p/\partial n < 0$ on $\partial\Omega$.

**Proof.** This a direct consequence of Theorem 3.5 in [8].

<figure><img src="image_12.png" /></figure>

In the discrete setting, the value of a harmonic function at a vertex is the mean value on the set of direct neighbors.

**Proposition 11.** Let $(V, E, c, \Gamma)$ be a resistive network, and $p$ a field that is harmonic on $\tilde{V}$. Then, for any $x \in \tilde{V}$, $p_x$ is an explicit convex combination of values at neighboring vertices:

$$p_x = \sum_{y \sim x} \pi_{xy} p_y, \text{ with } \pi_{xy} = \frac{c_{xy}}{C_x}, C_x = \sum_{y \sim x} c_{xy}.$$

**Proof.** This is a straight consequence of the harmonicity at $x$.

<figure><img src="image_13.png" /></figure>
Note that the counterpart of Proposition 7 for the continuous Laplacian is simply the definition of $L$, introduced as the difference between the vertex value and the mean value over neighbors.

**Remark 12.** The explicit formula given by the latter proposition is only valid for the strict neighborhood of a vertex. One can actually express a more general property, and express the value at some point $x$ as a convex combination of values taken in a set $\gamma$ of vertices that surrounds $x$ in some sense. Yet, unlike in the continuous setting, this combination is not fully explicit, it relies on the harmonic measure over $\gamma$, relative to $x$ (see Proposition 19).

### 2.3. Stochastic setting

The Dirichlet problem (4) has the following stochastic interpretations (see e.g. [7]).

**Proposition 13.** Consider a smooth and bounded domain $\Omega$ in $\mathbb{R}^d$, and $P$ a smooth function defined on the boundary $\Gamma = \partial\Omega$. For any $x \in \Omega$, let $w_x^t$ be the brownian motion from $x$. Let $\tau_x$ be the hitting time of $\Gamma = \partial\Omega$ by $w_x$, i.e.

$$ \tau_x = \inf\{t, w_x^t \in \Gamma.\} $$

Define $p$ as the expected value of $P$ at the hitting point:

$$ p(x) = \mathbb{E}(P(w_x^{\tau_x})). $$

Then $p$ is the solution to Dirichlet Problem (4).

**Proof.** See [7] (Chapter IX), or [9] (Chapter 4). $\square$

<figure><img src="image_14.png" /></figure>

Another property involves the so-called harmonic measure relative to a point.

**Proposition 14.** Consider a smooth and bounded domain $\Omega$ in $\mathbb{R}^d$, and $x \in \Omega$. Let $w_{x_0}(t)$ be the brownian motion from $x_0$ and $\tau_{x_0}$ the hitting time of $\Gamma = \partial\Omega$ by $w_{x_0}$ like in the previous proposition. Then $Y_{x_0} = w_{x_0}^{\tau_{x_0}} \in \Gamma$ is a random variable valued in $\Gamma$, the law of which is the measure $\mu_{x_0}$ on $\Gamma$. This measure admits a density with respect to the Lebesgue measure on $\Gamma$, which is $-\partial q/\partial n$, where $q$ is the solution to the homogeneous Dirichlet problem with singular right-hand side

$$ \left| \begin{array}{ll} -\Delta q = \delta_{x_0} & \text{in } \Omega \\ q = 0 & \text{on } \Gamma. \end{array} \right. \qquad (4) $$

In other words, for any measurable subset $A \subset \Gamma$

$$ \mathbb{P}(Y_x \in A) = - \int_A \frac{\partial q}{\partial n}. $$

The measure $\mu_{x_0}$ with density $-\partial q/\partial n$ is called the harmonic measure relative to $x_0$.

**Proof.** See [10]. $\square$

<figure><img src="image_15.png" /></figure>

As a direct consequence of the two previous propositions, one can express the solution to the Dirichlet problem

**Proposition 15.** Let $\Omega$ be a smooth and bounded domain, and $p$ the solution to Dirichlet problem (4) associated to the boundary condition $P$. Then, at any $x \in \Omega$, $p$ expresses

$$ p(x) = \int_{\partial\Omega} P(y) d\mu_x(y) $$

where $\mu_x$ is the harmonic measure relative to $x$.
This property is the conceptual core of Monte Carlo methods to solve non-homogeneous Dirichlet problems, we refer to the pioneering paper [11] for more details.

**Proposition 13** admits a discrete counterpart, where the brownian motion is replaced by the random walk canonically associated to the resistive network (see e.g. [12]).

**Proposition 16.** Let $(V, E, c, \Gamma)$ be a connected resistive network, with boundary $\Gamma \neq \emptyset$, and $P \in \mathbb{R}^\Gamma$ a collection of boundary values. We consider the Markov process on $V$ defined by $x \to y$ transition probabilities

$$ \pi_{xy} = \frac{c_{xy}}{C_x}, \quad C_x = \sum_{y \sim x} c_{xy}, $$

with $c_{xy} = 0$ whenever $(x, y) \notin E$. For any $x_0 \in V$, we denote by $(w_{x_0}^k)_k$ the random walk from $x_0$, according to the transition probabilities indicated above. We define the hitting time $k_{x_0}$ as

$$ k_{x_0} = \inf\{k, w_{x_0}^k \in \Gamma\} $$

and $Y_{x_0} = w_{x_0}^{k_{x_0}}$ the hitting point, that is a random variable valued in $\Gamma$. We finally define $p \in \mathbb{R}^V$ by

$$ p_{x_0} = \mathbb{E}(P_{Y_{x_0}}). $$

Then $p$ is the solution to the discrete Dirichlet problem (3).

**Proof.** If $x_0 \in \Gamma$, then $k_{x_0} = 0$ and $p_{x_0} = P_{x_0}$, so that $p$ verifies the boundary condition. If $x_0 \notin \Gamma$, the property is a straight consequence of the law of total expectation. Indeed,

$$ p_{x_0} = \sum_{y \sim x_0} \pi_{x_0 y} \mathbb{E}(P_{Y_y}) = \sum_{y \sim x_0} \pi_{x_0 y} p_y = \sum_{y \sim x_0} \frac{c_{x_0 y}}{C_{x_0}} p_y, $$

so that

$$ \sum_{y \sim x_0} c_{x_0 y} (p_{x_0} - p_y) = 0, $$

i.e. $p$ is harmonic, which ends the proof. $\square$

<figure><img src="image_16.png" /></figure>

**Proposition 14** also admits a discrete counterpart, which characterizes the discrete harmonic measure.

**Proposition 17.** (Harmonic measure in the discrete setting)

Let $(V, E, c, \Gamma)$ be a connected resistive network, with boundary $\Gamma \neq \emptyset$. Like in the previous proposition, for any $x_0 \in \mathring{V}$, we define $Y_{x_0} = w_{x_0}^{k_{x_0}}$ as the hitting point in $\Gamma$ of the random walk from $x_0$. For any $y \in \Gamma$, the probability that $Y_{x_0}$ is $y_0$ writes

$$ \mathbb{P}(Y_{x_0} = y_0) = \partial u(y_0) = -\partial c \partial^* p(y_0), $$

where $p$ is the solution to the Poisson problem

$$ \left| \begin{array}{ll} (\partial c \partial^* p)_x = \delta_{x_0, x} & \forall x \in \mathring{V} \\ p = 0 & \text{on } \Gamma. \end{array} \right. \qquad (5) $$

where $\delta_{x_0, x}$ is the Kronecker delta. The field $-\partial c \partial^* p$ on $\Gamma$ is the harmonic measure with respect to $x_0$.

**Proof.** The proof relies on an auxiliary Laplace problem with Dirichlet boundary conditions. Let $q$ be the solution to

$$ \left| \begin{array}{ll} \partial c \partial^* q(x) = 0 & \forall x \in \mathring{V} \\ q_y = \delta_{y, y_0} & \forall y \in \Gamma. \end{array} \right. \qquad (6) $$

Since $q$ is harmonic on $\mathring{V}$, and $p$ vanishes on $\Gamma$, it holds that

$$ \sum_{x \in \mathring{V}} p_x (\partial c \partial^* q)_x = 0. $$
Since the Laplacian $\partial c \partial^*$ is selfadjoint, the latter expression also writes

$$ \sum_{x \in V} p_x (\partial c \partial^* q)_x = \sum_{x \in V} q_x (\partial c \partial^* p)_x = q_{x_0} + (\partial c \partial^* p)_{y_0}, $$

so that $q_{x_0} = -(\partial c \partial^* p)_{y_0}$. Now, thanks to Proposition 16, $q_{x_0}$ is the probability that $Y_{x_0}$ is $y_0$, which ends the proof. □

**Remark 18.** In a euclidean domain, if one extends the pressure field of Prop. 14 by 0, $-\partial p/\partial n$ expresses mass unbalance on the boundary $\Gamma$. It is also $u \cdot n$, with $u = -\nabla p$. It identifies to the singular part of $-\Delta p$, which is a single layer distribution on $\Gamma$. In the discrete setting, there is no such notion as “normal derivative”, but the mass unbalance at a boundary vertex $y_0$ is $(\partial u)_{y_0}$ (that is the outflow through $y_0$), with $u = -c \partial^* p$, so that the mass unbalance writes $-(\partial c \partial^* p)_{y_0}$, which is therefore the discrete counterpart of $-\partial p/\partial n$ at $y_0$, as expressed by the previous proposition.

The notion of harmonic measure makes it possible to establish a generalization of the Mean Value Property in the discrete setting. The following proposition is a straight discrete counterpart of Prop. 15.

**Proposition 19.** Let $(V, E, c, \Gamma)$ be a connected resistive network, with boundary $\Gamma \neq \emptyset$, and $p \in \mathbb{R}^V$ a field that is harmonic in $\tilde{V}$. Let $x_0 \in \tilde{V}$ be given, and $\gamma \subset \tilde{V}$ a set of interior vertices which isolates $x_0$ from $\Gamma$, in the sense that all the path from $x_0$ to any vertex of $\Gamma$ necessarily meet $\gamma$. We denote by $V_{x_0}$ the set of vertices which can be connected to $x_0$ without meeting $\gamma$. Then $p_x$ writes as a convex combination of values taken by $p$ on $\gamma$:

$$ p_x = \sum_{y \in \gamma} \theta_y p_y, \quad \theta_y \ge 0 \quad \forall y \in \gamma, \quad \sum_{y \in \gamma} \theta_y = 1. \qquad (7) $$

The field of weights $(\theta_y)_\gamma$ is the harmonic measure on $\gamma$ relative to $x_0$, in the sense of Proposition 17, i.e.

$$ \theta_y = -(\partial c \partial^* q)_y, $$

where $q$ is the solution to the Dirichlet problem on $V_{x_0} \cup \gamma$

$$ \left| \begin{array}{l} (\partial c \partial^* q)_x = \delta_{x_0, x} \quad \forall x \in \tilde{V}_{x_0} \\ q_y = p_y \quad \forall y \in \gamma. \end{array} \right. \qquad (8) $$

**Proof.** Let us consider the subnetwork $(V', E', c', \Gamma')$, where $V' = V_{x_0} \cup \gamma$, $E'$ contains the edges of $E$ that connect two vertices in $V'$, $c'$ is simply the restriction of $c$ to $E'$, and $\Gamma' = \gamma$. We apply Proposition 17 to the subnetwork $(V', E', \Gamma')$. It ensures that $\theta_y = -(\partial c \partial^* q)_y$ is the probability that a random walk from $x_0$ meets $\gamma$ at $y$, where $q$ solves (8). Now since $\gamma$ isolates $x_0$ from $\Gamma$, any random walk starting from $x_0$ and ending at first hit with $\gamma$ does not meet $\Gamma$, thus it only visits vertices in $\tilde{V}$, at which the field $p$ is harmonic, which makes it possible to use Proposition 16 applied to the new network $(V', E', \Gamma')$, with $\Gamma' = \gamma$. As a consequence, $p_{x_0}$ is the expected value of $p$ at the first hit with $\gamma$, which is expressed by (7). □

### 3. Macroscopic and microscopic congestion models

The approach which we present in this section was initially motivated by the handling of collisions in fluid-particle simulations (see e.g. [13–15]). We shall present a common strategy to handle congestion at microscopic and macroscopic scales. This strategy can be carried out for second order in time evolution problems, see [15] for microscopic granular flows, and [16] for the macroscopic version. We also refer to [17], where the authors propose a numerical scheme to solve Euler equations with a singular pressure meant to account for congestion, or [18] for the asymptotic analysis of a similar system. We shall present here the overdamped versions of these
inertial problems, with a hard handling of the maximal density constraint, in the context of crowd motion modeling.

We describe in this section the macroscopic and the microscopic settings, and we show how they both lead to a Poisson-like problem. The next section proposes a deeper look on the discrepancy between both settings, and a detailed analysis of the pathologies of the discrete Laplacian.

### Macroscopic setting

The macroscopic crowd motion model proposed in [19] reads as follows. Given a domain $\mathcal{O}$ (which covers the zone possibly occupied by the crowd), the crowd is represented by a density $\rho = \rho(x, t)$ subject to remain below a certain maximal value, which we set to 1. The set of feasible density is thus

$$K = \{\rho \in L^1(\mathcal{O}), 0 \le \rho \le 1 \text{ a.e.}\}.$$

We denote by $U = U(x)$ the desired velocity, that is given: $U(x)$ is the velocity that a person at $x$ would take, would they be alone. We write that $\rho$ is transported by the actual velocity $u$:

$$\partial_t \rho + \nabla \cdot (\rho u) = 0, \qquad (9)$$

where $u$ is the closest to $U$ among all those velocities which do not lead to a violation of the constraint, i.e. which are expansive in the zone which are already saturated. More formally, it amounts to prescribe a nonnegative divergence in the saturated zone. This constraint can be expressed in a dual way by subjecting $u$ to belong to

$$C_{\rho} = \left\{ v \in L^2(\mathcal{O}), \int_{\mathcal{O}} v \cdot \nabla q \le 0 \quad \forall q \in H_{\rho}^1(\mathcal{O}), q \ge 0 \text{ a.e.} \right\}$$

were

$$H_{\rho}^{1}(\mathcal{O}) = \{q \in H^{1}(\mathcal{O}), q(1-\rho) = 0 \text{ a.e.}\}.$$

Note that the condition $q(1 - \rho) = 0$ contrains $q$ to vanish in zones where $\rho < 1$, i.e. where the constraint is not saturated. We may now write the full model by defining the actual velocity as

$$u = P_{C_{\rho}}(U), \qquad (10)$$

where the projection is defined in the $L^2$ sense.

**Remark 20.** Model (10), like its microscopic counterpart (14) below, is not meant to account for sophisticated individual tendencies, it merely encodes a very basic strategy in terms of behavior: individual strive to achieve their desired velocity, regardless of others. The resulting interactions are physical in nature (contact forces), and they may induce nonlocal effects over connected saturated zones.

Boundary conditions can be accounted for through the definition of the pressure field. As an archetypal example, consider the situation represented in Figure 4 where $\mathcal{O}$ is a room with boundary $\Gamma$, union of $\Gamma_w$ (walls) and $\Gamma_{out}$ (exit). The saturated zone $[\rho = 1]$ corresponds to the grey domain, which we denote by $\Omega$. In this situation, we prescribe $q = 0$ on $\Gamma_{out}$, whereas $q$ is “free” on $\Gamma_w$. Note that, due to the constraint $q(1 - \rho) = 0$, the pressure vanishes on the upstream side of the saturated zone $\Omega$. Let us now comment on the definition of $C_\rho$, in this particular situation. Suppose $\Omega = [\rho = 1]$ is a smooth domain, and suppose $u \in C_\rho$ is smooth, considering smooth and nonnegative test functions $q$ in $H_\rho^1$, compactly supported in $\Omega$, the dual condition yields

$$\int_{\mathcal{O}} u \cdot \nabla q = - \int_{\Omega} q \nabla \cdot u \le 0,$$
<figure><img src="image_17.png" /><figcaption>**Figure 4.** Evacuation from a room.</figcaption></figure>

so that $\nabla \cdot u \ge 0$ a.e. in $\Omega$. Now consider tests functions which are no longer compactly supported in $\Omega$, which can take positive values on $\Gamma_w$. The same integration by part yields a boundary term, so that we obtain

$$\int_{\Gamma_w} q u \cdot n \le 0,$$

which imposes $u \cdot n \le 0$ on $\Gamma_w$: people may not exit through the wall.

We refer to [19] for a full analysis of the evolution problem (9)(10). This analysis is delicate due to the nonlinearity and lack of smoothness of $\rho \mapsto u = P_{C_\rho} U$, which rules out standard PDE approaches. The natural framework is the Wasserstein space, built upon a distance between densities that is based on optimal transportation. We shall disregard here these aspects, and rather focus on the projection step (10), which defines $u$. In a paradoxical way, we shall see that this step is much more delicate and rich at the microscopic scale than at the macroscopic one. Indeed, implementing (10) amounts to minimize the functional

$$v \mapsto J(v) = \int_{\mathcal{O}} |v - U|^2,$$

on the feasible set $C_\rho$. Since $C_\rho$ is closed and convex in $L^2(\mathcal{O})$, and since $J$ is continuous and strongly convex, the projection problem admits a unique solution $u$.

*Saddle-point formulation and Poisson problem*

Let us introduce

$$B: v \in L^2(\mathcal{O}) \longrightarrow -\nabla \cdot v \in H_{\rho}^{-1} = (H_{\rho}^{1})'$$

defined by

$$\langle -\nabla \cdot v, q \rangle = \int_{\mathcal{O}} v \cdot \nabla q.$$

The projection problem (10) can be formulated as a saddle-point problem (see e.g. [20]) for the Lagrangian

$$L(v, q) = J(v) + \langle Bv | q \rangle = \int_{\mathcal{O}} |v - U|^2 + \int_{\mathcal{O}} v \cdot \nabla q. \quad (11)$$

**Proposition 21.** *The Lagrangian defined by (11) admits a unique saddle point $(u, p)$ in $V \times \Lambda_+$, with*

$$V = L^2(\mathcal{O}), \Lambda_+ = \{q \in H_{\rho}^1(\mathcal{O}), q \ge 0 \text{ a.e. .}\}$$
This saddle point is characterized by $u \in L^2(\mathcal{O})$, $p \in H_p^1(\mathcal{O})$

$$ 
\left\| 
\begin{array}{ll}
u + \nabla p & = U \text{ in } \mathcal{O} \\
-\nabla \cdot u & \le 0 \text{ on } [\rho = 1] \\
p & \ge 0 \text{ a.e. in } \mathcal{O} \\
\int_{\mathcal{O}} u \cdot \nabla p & = 0
\end{array}
\right. 
\qquad (12) $$

where the second equation is meant in a weak sense, i.e.

$$ \int_{\mathcal{O}} u \cdot \nabla q \le 0 $$

for any non-negative test function which vanishes outside $[\rho = 1]$. The primal part of this saddle-point is the projection of $U$ on $C_p$.

**Proof.** This problem fits into the abstract framework of Prop. 32 in Sect. 5, in the infinite dimensional situation. Well-posedness comes from that fact that $B$ is surjective, which is a straight consequence of Poincaré inequality (see Proposition 3 in [21] for details). $\square$

Now consider a situation where the saturated zone is a smooth domain (like in the evacuation situation represented in Figure 4). As detailed in [21], if the desired velocity field in concentrating, i.e. $\nabla \cdot U < 0$, then $p$ cannot vanish on an open subset on $\Omega$ (it would rule out the constraint in this subdomain). Thus, the complementarity slackness condition (4-th line in (12)) imposes that the constraint is saturated, i.e. $\nabla \cdot u = 0$, in $\Omega$. As a consequence, one may eliminate the velocity to obtain a Poisson problem on the pressure with mixed boundary conditions:

$$ 
\left| 
\begin{array}{ll}
-\Delta p = -\nabla \cdot U & \text{in } \Omega \\
p = 0 & \text{on } \Gamma_{\text{out}} \cup \Gamma_{\text{up}} \\
\frac{\partial p}{\partial n} = 0 & \text{on } \Gamma_w.
\end{array}
\right. 
\qquad (13) $$

As we shall see in the next section, this Poisson problem is central in the behavior of this congestion model.

### Microscopic setting

The microscopic counterpart of the previous approach was introduced in [22]. It is Lagrangian in nature, based on the position of $N$ individuals $x = (x_1, \dots, x_N) \in \mathbb{R}^{2N}$. The individuals are identified to rigid discs of radius $r$, they are subject to a non-overlapping constraint, i.e. $g$ must remain in the feasible set

$$ K = \{x \in \mathbb{R}^{2N}, D_{ij}(x) = |x_j - x_i| - 2r \ge 0 \quad \forall i \neq j .\} $$

We consider a collection of desired velocities $U = (U_1, \dots, U_N) \in \mathbb{R}^{2N}$. The constraint on the divergence at the macroscopic level becomes here a non-overlapping constraint: Feasible velocities are those fields which do not decrease distances that are already 0. For any $x \in K$, the feasible set for velocities writes

$$ C_x = \{v \in \mathbb{R}^{2N}, D_{ij}(x) = 0 \implies e_{ij} \cdot (v_i - v_j) \le 0\}, $$

where $e_{ij} = (x_j - x_i) / |x_j - x_i|$ is the unit vector pointing from $x_i$ to $x_j$. The evolution problem, which is the microscopic counterpart of (9)(10), writes

$$ 
\left| 
\begin{array}{l}
\displaystyle \frac{dx}{dt} = u \\
u = P_{C_x} U.
\end{array}
\right. 
\qquad (14) $$
<figure><img src="image_18.png" /><figcaption>**Figure 5.** One-dimensional cluster</figcaption></figure>

Note that, in this purely Lagrangian setting, the Eulerian transport equation (9) is replaced here by a simple definition of the velocity as the time derivative of the position vector.

Like in the macroscopic context, the set of feasible velocities is a closed convex cone, which ensures existence and uniqueness of a projection. The analysis of the time evolution problem, which is not the topic of the present paper, is detailed in [22], it relies on the so-called prox-regularity of the position feasible set $K$, which essentially says that the projection on this set is locally well-defined. We shall focus here on the projection step which defines $u$. Let us first introduce a matrix formulation of the constraints : the set of feasible velocities can be written

$$C_x = \{v, Bv \le 0\}$$

where each row of $B$ stands for an active contact, i.e. a couple $(i, j)$ such that $D_{ij} = 0$, and thus writes

$$B_{ij} = (0, 0, \dots, e_{ij}, 0, \dots, 0, -e_{ij}, 0, \dots, 0) \in \mathbb{R}^{2N}.$$

The matrix $B$ belongs to $\mathcal{M}_{N_c \times 2N}(\mathbb{R})$, where $N_c$ is the current number of contacts. This number may vary in time, but we shall here focus on the instantaneous problem.

The constrained minimisation problem admits a saddle-point formulation for the Lagrangian

$$L(v, q) = J(v) + \langle Bv | q \rangle = \frac{1}{2} \sum_{i=1}^{N} |v_i - U_i|^2 + \sum_{i<j, D_{ij}=0} p_{ij} \times e_{ij} \cdot (v_i - v_j), \quad (15)$$

as expressed in the following proposition (discrete counterpart of Prop. 21).

**Proposition 22.** *The Lagrangian defined by (15) admits a saddle point $(u, p)$ in $V \times \Lambda_+$, with*

$$V = \mathbb{R}^{2N}, \Lambda_+ = \mathbb{R}_+^{N_c}.$$

*This saddle point is characterized by*

$$\left\{
\begin{array}{l}
u + B^* p = U \\
Bu \le 0 \\
p \ge 0 \\
\langle p | Bu \rangle = 0.
\end{array}
\right.
\qquad (16)$$

*The primal part of this saddle-point is the projection of $U$ on $C_x$.*

**Proof.** This elementary property relies on classical arguments: like in the macroscopic setting (Prop. 21), this problem fits into the abstract framework of Prop. 32, in the case where the number of constraints is finite. Note that no assumption is made on $B$, in particular surjectivity is not required, so that uniqueness may be ruled out for the dual component $p$. □
All matrices take a simple form in the one-dimensional setting. Consider a cluster of $N$ discs in a row (see Fig. 5). Matrix $B$, discrete opposite of the divergence, and $BB^*$, write

$$ B = \begin{pmatrix} 1 & -1 & 0 & \dots \\ 0 & 1 & -1 & \dots \\ 0 & 0 & \ddots & \dots \\ 0 & 0 & \dots & 1 & -1 \end{pmatrix}, \quad BB^* = \begin{pmatrix} 2 & -1 & 0 & \dots & 0 \\ -1 & 2 & -1 & 0 & \dots \\ \vdots & \vdots & \ddots & \ddots & \vdots \\ -1 & 2 & -1 & 0 & \dots \\ 0 & \dots & 0 & -1 & 2 \end{pmatrix}, \qquad (17) $$

which is the discrete Laplacian in dimension one.

**Remark 23.** The operator $B$ plays here the role of $\partial$ for resistive networks (see Section 2.1 and Fig. 3), which is the discrete counterpart of the opposite of the divergence. Mass balance is replaced here by overlap estimate. Indeed, by definition, each line of $B$ expresses the first order term of the expansion of the distance, e.g

$$ D_{ij}(x + \epsilon U) = \underbrace{D_{ij}(x)}_{=0} + \epsilon(BU)_{ij} + o(\epsilon). $$

Accordingly, $B^*$ plays the role of $\partial^*$, which is a discrete gradient. For a resistive fluid network, a pressure field induces a collection of fluxes proportional to $\partial^* p$, which write $u_{xy} = c_{xy}(p_x - p_y)$ in each edge $(x, y)$. In the present context, a pressure field $p \in \mathbb{R}^{N_c}$ induces a force field $B^* p$ over the primal degree of freedom, which are the positions of the grains. Note that the fact that $B$ and $B^*$ are mutually adjoint is not fortuitous like in the case of resistive networks, this property is inherent to the saddle-point formulation which defines $p$.

**Remark 24.** Note that this matrix corresponds to the Finite Difference discretization of a Laplace problem with Dirichlet boundary condition. Since the grain 1 (see Fig. 5) has no neighbor on its left-hand side, it is as if a zero pressure were imposed at the grey circle on its left. In the Finite Difference context, these grey circles correspond to boundary points at which a zero value is prescribed, so that they do not correspond to unknowns, and they are in general removed from the matrix formulation. This will play an important role in what follows: we shall define Laplace-like operators on clusters of discs, and those operators will act on pressures which are defined on edges of the contact network, which is the dual graph. The resulting discrete Laplacian corresponds to a situation where a zero pressure is prescribed on the border.

For discs in $\mathbb{R}^2$, following the idea that $B$ is the discrete counterpart of the opposite of the divergence, it can be expected that the matrix $BB^*$ can be interpreted as a discrete Laplacian. Figure 6 represents such a cluster of discs, the associated primal graph, and the associated dual graph (in red), the vertices of which are the points of application of the interactions forces. The operator $BB^*$ acts on the dual graph of the contact network. As illustrated in the figure, the dual graph is not planar as soon as some discs have more than 3 neighbors, whereas the primal graph is planar. The thick red lines in the figure correspond to the stencil associated to a “vertex” of the dual graph, that is the contact $(i, k)$. This dual vertex has 7 neighbors, which are the other contacts that implicate either $i$ or $k$. One may imagine that some substance circulates among discs, with a flux between two discs given by a Fick-like law. This would lead to a standard Laplacian of the type $\partial_c \partial^*$, defined on the primal graph (in black in the figure), like in Section 1. The maximal degree would be 5 (grain $k$ has 5 neighbors), whereas the maximal degree of the dual graph is 7 here. As we shall see, the operator $BB^*$ which emerges from the saddle-point formulation, which is dual in some way to $\partial_c \partial^*$, presents very different characteristics.

Like in the macroscopic setting, one may deduce from Prop. 16 that the pressure $p$ verifies some sort of Poisson problem. Yet, an extra difficulty comes from the fact that the constraint
<figure><img src="image_19.png" /><figcaption>**Figure 6.** Primal and dual graphs</figcaption></figure>

$Bu \le 0$ is not necessarily saturated, even in the case of a concentrating velocity field $U$, i.e. such that $BU > 0$, which is the first major difference with the macroscopic setting. Let us explain this first discrepancy in a very simple situation, illustrated by Figure 7 (right), which is a discrete version of the macroscopic situation represented in the left-hand side of the figure. The situation is represented in the upper right corner of the figure : 4 discs form a cluster, and the desired velocities push them toward the center. For any two discs in contacts, moving the discs along desired velocities would lead to a violation of the constraint, which writes $BU > 0$.

A macroscopic counterpart of this situation is represented on the left-hand part of the figure: the domain is an ellipse, and the desired velocity points toward the center, in such a way that all parts of the domain are compressed by $U$. It writes $-\nabla \cdot U > 0$. In the macroscopic situation, the pressure field is the solution to the Poisson problem $-\Delta p = -\nabla \cdot U > 0$, so that $p > 0$ in $\Omega$ by Corollary 10. In this situation, minimizing $|v - U|$ under the constraint $-\nabla \cdot v \le 0$ thus leads to a pressure field that is positive, so that minimizing over divergence-free fields would lead to the very same solution. In the microscopic setting the situation is fully different. Minimizing over $\text{ker } B$ amounts to consider that grains are glued to each other (they actually can slip on each other, but distances of grains in contact have to remain equal to 0). For symmetry reasons, the projection of $U$ over $\text{ker } B$ is obviously 0. In the saddle-point formulation, the interaction force between the two grains that are on the vertical axis is *negative*, whereas the other pressures are positive (see the '+' and '-' signs on the figure), which is in full opposition to the macroscopic situation. As we shall see, it comes from the fact the Laplace operator that we are going to exhibit does not verify any maximum principle. If one now considers the inequality constraint $Bu \le 0$, the solution will be different: the two grains on the vertical axis shall be pushed apart, due to the horizontal squeezing action on the grains on the horizontal axis (*wedge* effect).

The example that we just considered shows that supposing $BU > 0$ (grains are pushed against each other) does not guarantee that the constraint is saturated (i.e. $Bu = 0$). To obtain a Laplace operator by eliminating the velocity thus requires that we first remove from $B$ the rows which correspond to non-active constraints, i.e. such that the actual velocity fields create a *positive* relative velocity: grains are strictly pulled apart each other.
<figure><img src="image_20.png" /><figcaption>**Figure 7.** Wedge principle.</figcaption></figure>

**Proposition 25.** For a given desired velocity field $U$, $(u, p)$ being a solution to the saddle-point formulation (16), we denote by $\mathring{B}$ the matrix obtained from $B$ by removing all the rows associated to a non active constraint, i.e. for couple $(i, j)$ such that $e_{ij} \cdot (v_i - v_j) < 0$. Then $p$ verifies the discrete Poisson problem

$$ \mathring{B}\mathring{B}^*\mathring{p} = \mathring{B}U, \qquad (18) $$

where $\mathring{p}$ is obtained from $p$ by removing all the entries that correspond to removed constraints.

**Proof.** By construction, we have that $\mathring{B}u = 0$ (all the rows with a strict inequality have been removed). Now, thanks to the complementarity condition, for all those rows $(i, j)$ which have been removed, it holds that $p_{ij} = 0$, so that $\mathring{B}^*\mathring{p} = B^*p$. We thus obtain the discrete Darcy problem

$$ \left| \begin{array}{l} u + \mathring{B}^* \mathring{p} = U \\ \mathring{B} u = 0 \end{array} \right. \qquad (19) $$

which yields, after elimination of the velocity, to the discrete Poisson equation (18). $\square$

Let us insist on the fact that the discrete Laplacian $\mathring{B}\mathring{B}^*$ does not only depend on the configuration. For a given $x$ in $K$ (collection of non-overlapping discs), $BB^*$ is the basic discrete Laplacian but, given a compressing field $U$ (i.e. such that $BU > 0$), the induced pressure will verify a Poisson problem $\mathring{B}\mathring{B}^*\mathring{p} = \mathring{B}U$, where $\mathring{B}$ is obtained by removing some rows from $B$, in a way that depends on $U$. In other words, there is no canonical Laplacian associated to the dual graph, but a (finite) collection of Laplacian operators which are defined on subgraphs of the dual graph.

**Remark 26.** Let us give a very concrete illustration of the previous considerations. Consider a collection of many coins on a table, with no overlapping but possibly touching. If one exerts centripetal forces from the outside of the cluster, the induced motion may push some coins away from each other, and the collection of couples that are split depends on the way the squeezing force is prescribed.
<figure><img src="image_21.png" /><figcaption>**Figure 8.** Acute angle condition</figcaption></figure>

To better understand the origin of this pathology, let us describe in more detail the matrix $B$, in the form of a proposition which expresses a negative property.

**Proposition 27.** Let $x \in K$ represent a collection of $N$ non-overlapping discs. If there exists at least a disc $i$ which is in contact with $j$ and $k$ such that $e_{ij}$ and $e_{jk}$ form an acute angle, then the matrix $BB^*$ is not a $Z$-matrix, i.e. it contains extra-diagonal positive entries.

**Proof.** This comes from the expression of the matrix $BB^*$. The simple example of such a situation is represented in Fig. 8, which is an extraction of a 3-cluster from the configuration represented in Fig. 6. Let us compute the $2 \times 2$ contribution to $BB^*$ that comes from the interaction between $(i, j)$ and $(i, k)$ (red dots in the figure). The corresponding matrices are

$$B_{(i,j),(i,k)} = \begin{pmatrix} e_{ij} & -e_{ij} & 0 \\ e_{ik} & 0 & -e_{ik} \end{pmatrix} \in \mathcal{M}_{2,6}(\mathbb{R}),$$

so that the contribution of this interaction to the global Laplacian is

$$B_{(i,j),(i,k)} B_{(i,j),(i,k)}^* = \begin{pmatrix} 2|e_{ij}| & e_{ij} \cdot e_{ik} \\ e_{ij} \cdot e_{ik} & 2|e_{ik}| \end{pmatrix} = \begin{pmatrix} 2 & 1/2 \\ 1/2 & 2 \end{pmatrix} \in \mathcal{M}_{2,2}(\mathbb{R}).$$

By gathering all those contributions, the extra diagonal term which corresponds to the $(i, j)$, $(i, k)$ interaction will remain unchanged, thus positive. As a consequence, the resulting matrix $BB^*$ has positive extra-diagonal entries. $\square$

In the continuous setting, a smooth function defined in a smooth connected domain is constant if and only if its gradient is identically 0. In the discrete setting, where $B^*$ plays the role of the gradient, this equivalence is no longer true, and this negative property is not due to boundary effects.

If one considers the situation represented in Fig. 6, the disc $k$ is in contact with 5 discs, and the contact points are not distributed in a symmetric way. As a consequence, a unit pressure field induces a non-zero force on grain $k$, i.e. $B^*1 \neq 0$.

Conversely, in highly congested situations, the discrete gradient $B^*$ may have a rich kernel. Consider for example a triangular lattice of discs with a common radius. If the number of discs is $N$, there are $2N$ primal degrees of freedom, whereas the number of contacts $N_c$ scales like $3N$.
The kernel of $B \in \mathcal{M}_{N_c, 2N}$ has therefore a dimension that scales like $N$. As a direct consequence, the Poisson problem

$$BB^* p = b$$

which is the discrete counterpart of a standard Poisson problem with Dirichlet boundary conditions (see Remark 24), is highly *ill-posed* in terms of uniqueness. Let us nevertheless add that, even in situations where the problem is very degenerate, the unilateral homogenous problem is well-posed: the problem $BB^* p = 0$, with $p \in \mathbb{R}_+^{N_c}$ admits 0 as unique solution. This is a straight consequence of Hahn-Banach theorem, as detailed in [15]. A simple illustration of this property is the following: consider a static collection of coins on a table, possibly in contact with each other. Assuming that coins are not glued to each other, the question is: is it possible to determine interaction forces ? The answer is yes. Consider an extreme coin (the center of which is not in the convex hull of other centers). Since forces are repulsive, if the forces acting upon this coin were not zero, they would push the coin away from the cluster. So they are zero, the coin can be removed without perturbing the equilibrium, and the process can be further implemented until exhaustion. Note that, in the same situation, if some coins are glued together (i.e. pressures can be positive or negative), it is impossible to determine interaction forces from the sole equilibrium property. It can also be proved that the solution set of the Poisson problem with positivity constraint is bounded. Indeed (see again [15], Lemma 1), it holds that $ker B^* \cap \mathbb{R}_+^{N_c}$ is bounded, so that any Poisson problem $BB^* p = b$ set in $ker B^* \cap \mathbb{R}_+^{N_c}$ admits a bounded (possibly empty) solution set.

Before exploring the effect of the (negative) properties of $BB^*$ upon the behavior of the evolution model (14), let us summarize these properties. Its form provides it with a variational setting. Indeed, as illustrated in Fig. 14, the pressure that is a solution to the saddle-point formulation (16), minimizes $|U - B^*q|^2$. As for the other properties listed in Section 2, most of them are *not* verified by $BB^*$. Since some extra diagonal terms can be positive, there is no mean value properties (Prop. 11 is not true), no stochastic interpretation (Props. 16, 17, and 19 do not make any sense, since there is no random walk canonically associated to $BB^*$).

## 4. Capacity drop and Faster is Slower effect at the macroscopic and macroscopic scales

We explain in this section how the “nice” properties of the macroscopic (standard) Laplace operator prevent the corresponding model to reproduce some observed effect and, on the other side, how the pathologies of the microscopic operator make it possible to reproduce these effects.

### Capacity drop

In the context of evacuation, the *Capacity Drop* (CD) effect designs a reduction of the flux though an exit door or a bottleneck when the number of evacuees upstream increases, see e.g. [23]. Let us first show that the macroscopic model (9) (10) does *not* reproduce the CD effect.

**Proposition 28.** Let us consider the situation represented in Fig. 4: pedestrians are represented at some time by a density $\rho$ upstream the exit $\Gamma_{out}$, the saturated zone ($\rho = 1$) is a smooth domain $\Omega$, and the desired velocity field points to the door, in such a way that $\nabla \cdot U \le 0$, and this quantity is not identically 0 on $\Omega$. Then, according to the model (9) (10), people actually exit faster than they would if there were no congestion.

**Proof.** This is a straight consequence of the Maximum Principle. We give here a formal proof of this property, and we refer to [21] for a more detailed proof. Since $U$ is concentrating, the constraint $\nabla \cdot u \ge 0$ is saturated over $\Omega$, so that $\nabla \cdot u$ is identically 0. The unilateral Darcy problem (12) thus becomes a standard Darcy problem. Elimination of the velocity yields the Poisson problem

$$-\Delta p = -\nabla \cdot U,$$
<figure><img src="image_22.png" /><figcaption>**Figure 9.** Principle of a static jam</figcaption></figure>

with homogeneous Dirichlet boundary conditions on $\Gamma_{out}$ (the exit), $\Gamma_{up}$ (the upstream front), and homogeneous Neumann condition on $\Gamma_{w}$ (walls). The field $p$ is such that $-\Delta p \ge 0$, and not identically 0 in $\omega$, so that $p > 0$ in $\Omega$. As a consequence, for any $x \in \Gamma_{out}$ the value of $p$ is strictly less that $p(y)$ for any $y \in \Omega$ and, by Corollary 10, $\partial p/\partial n(x) < 0$. As a consequence

$$u \cdot n = (U - \nabla p) \cdot n = U \cdot n - \frac{\partial p}{\partial n}(x) > U \cdot n,$$

which exactly states that people exit faster than they would if they were alone (i.e. if they adopted their desired velocity). □

The previous property shows that the macroscopic model in its basic form is not relevant to recover the Capacity Drop phenomenon. Note that some attempts have been proposed to incorporate this phenomenon in macroscopic models. In [24, 25], the authors prescribe a limitation of the flux at the exit, in a way that non-locally depends on the number of people accumulated upstream the exit, which makes it possible to enforce the CD phenomenon. Yet, up to our knowledge, no macroscopic model is currently able to natively reproduce this phenomenon, in a way that would explain it.

We shall investigate now how the microscopic setting, because of the very defects of the discrete Laplacian $BB^*$, reproduce the Capacity Drop phenomenon. The discrete model actually exhibits an extreme version of the Capacity Drop phenomenon, that is the spontaneous emergence of static jams.

**Definition 29.** (Static jam)

In an evacuation situation like the one represented in Fig. 4, a static jam is an equilibrium point of the evolution problem (14), i.e. a configuration $x \in K$ such that $P_{C_x} U = 0$. By Proposition 16, it can also be defined as any configuration $x \in K$ such that there exists a pressure field $p \in \mathbb{R}_+^{N_c}$ which balances $U$, in the sense that $B^* p = U$.

We may express in an informal way the ability of the microscopic model to reproduce the Capacity Drop phenomenon:
<figure><img src="image_23.png" /><figcaption>**Figure 10.** Stable static jam</figcaption></figure>

In an evacuation situation like the one represented in Fig. 4, with desired velocities pointing to the outside of the room through the exit door, the microscopic evolution model (14) leads in some situations to static jams in the sense of Definition 29, even when the width of the exit is significantly larger than the size of an individual.

The mechanism of jams is illustrated in Fig. 9. The bottom corresponds to the macroscopic situation, it illustrates Proposition 28: the red arrows correspond to the variations of velocities due to congestion, i.e. $-\nabla p$, which points outward the domain, meaning that people exit faster than they would if they were alone (capacity rise phenomenon). This property is a direct consequence of the Maximum Principle for $-\Delta$, as detailed in Prop. 28. The top part of the figure corresponds to the microscopic situation. According to Remark 24, one can consider that the grey dots, located at the zone of grains facing the exit, correspond to points where a zero pressure is prescribed. Now consider positive pressures between grain in contact, the corresponding correction on the velocity, that is $-B^*p$ (discrete counterpart of $-\nabla p$), points inward the domain, meaning that people exit slower than they would if they were alone (CD phenomenon). The previous example explains the mechanism through which the desired velocity can be balanced by interaction forces. In can be checked that such jams spontaneously appear in numerical simulations, as illustrated by Fig. 10, which represents a jammed situation with the desired velocities in black, and the forces induced by pressure in blue.

Figure 11 represents the pressure field during an evacuation process. Pressures are associated to edges between two centers, and the width of edges are proportional to pressure values. Pressure arches are quite visible, in particular at the very exit, which instantiates the general principle described in Fig. 9. Note also that the highest pressures implicate discs upfront, to be compared to the macroscopic situation, where maximal pressures are attained on a zone significantly upstream the exit, so that people at the exit are pushed ahead.

### Faster is Slower effect

The Faster is Slower (FiS) effect pertains to a wide class of systems which may perform globally worse as their individual components strive to do better (see e.g. [26]). In the context of
<figure><img src="image_24.png" /><figcaption>**Figure 11.** Macroscopic (top) and microscopic (bottom) pressure fields.</figcaption></figure>

evacuation, it corresponds to a situation where the eagerness of some individuals to exit a room leads to a decrease of the evacuation efficiency. We explain here why the macroscopic model (9) (10) does *not* make it possible to recover the FiS effect, whereas its microscopic counterpart (14) reproduces this effect. We propose a partial explanation of this paradoxical phenomenon, which relies on the *defects* of the underlying Laplacian.

We consider again the situation represented by Fig. 4, with a desired velocity field $U$ which points toward the exit. In [27] it was proposed to change $U$ by introducing a scalar speed
correction field $\beta(x)$, $\beta > 1$ corresponding to a speed increase, and to investigate the dependence of the flux through the exit upon $\beta$. We extend here this approach by allowing variations of $U$ in any direction, and we shall keep $\beta$ to denote the variation around $U$, so that $U + \beta$ is the modified velocity field, and we shall investigate the dependence of the flux upon $\beta$. The FiS issue formalizes in this context:

Are there any zones such that a variation of $U$ in its own direction (so that $\beta \cdot U > 0$) induces a *decrease* of the flux at the exit?

The pressure $p_\beta$ associated to this desired velocity $U + \beta$ is the solution to the Poisson problem

$$ -\Delta p_{\beta} = -\nabla \cdot U - \nabla \cdot \beta, \quad (20) $$

with homogeneous Dirichlet boundary conditions $p_\beta = 0$ on $\Gamma_{out}$ and $\Gamma_{up}$, and homogeneous Neuman boundary conditions $\partial p_\beta / \partial n = 0$ on $\Gamma_w$. The associated flux through the exit is thus

$$ J(\beta) = \int_{\Gamma_{out}} u_{\beta} \cdot n = \int_{\Gamma_{out}} U \cdot n - \int_{\Gamma_{out}} \frac{\partial p_{\beta}}{\partial n}. \quad (21) $$

**Proposition 30.** Let us assume that the desired velocity $U$ is tangent to the wall. Let $\beta \mapsto J(\beta)$ the flux functional defined by (21). We assume that desired velocities are not modified at the exit, i.e. $\beta = 0$ on $\Gamma_{out}$. The gradient of $J$ at $\beta$, that is the field $g$ such that

$$ J(\beta + \tilde{\beta}) = J(\beta) + \int_{\Omega} g \cdot \tilde{\beta} + o(\tilde{\beta}), $$

for all variations $\tilde{\beta}$ vanishing on $\Gamma_{out}$ is $-\nabla q$, where $q$ is the solution to the adjoint problem

$$ \left| \begin{array}{ll} -\Delta q = 0 & \text{in } \Omega, \\ q = -1 & \text{on } \Gamma_{out} \\ q = 0 & \text{on } \Gamma_{up} \\ \\ \frac{\partial q}{\partial n} = 0 & \text{on } \Gamma_w. \end{array} \right. \quad (22) $$

**Proof.** The gradient is determined by the adjoint problem method. We introduce a dual variable $q$, that is a function defined in $\Omega$, and which makes it possible to express the relation between the control variable $\beta$ and the state variable $p$. For any such $q$ (assumed to be smooth), from (20), we have that

$$ b(p, q) = \int_{\Omega} \nabla p \cdot \nabla q - \int_{\Gamma} \frac{\partial p}{\partial n} q + \int_{\Omega} q \nabla \cdot U + \int_{\Omega} q \nabla \cdot \beta = 0, \quad (23) $$

where the boundary integral over $\Gamma$ can be replaced by an integral over $\Gamma_0 = \Gamma_{out} \cup \Gamma_{up}$, thanks to the boundary condition on $\Gamma_w$. Now, since we assume that desired velocities are not modified at the exit, and since we are interested in variations of the flux only, we may drop the first term of (21) in the definition of $J$:

$$ J(\beta) = - \int_{\Gamma_{out}} \frac{\partial p_{\beta}}{\partial n}. $$

The Lagrangian

$$ L(p, \beta, q) = - \underbrace{\int_{\Gamma_{out}} \frac{\partial p}{\partial n}}_{\text{flux}} + \underbrace{\int_{\Omega} \nabla p \cdot \nabla q - \int_{\Gamma_0} \frac{\partial p}{\partial n} q + \int_{\Omega} q \nabla \cdot U + \int_{\Omega} q \nabla \cdot \beta}_{b(p,q)} $$

is such that, for any $p_\beta$ coupled with $\beta$ by the state equation, $b(p_\beta, q) = 0$ for any $q$. Thus, for any $\beta$ and any $q$,

$$ L(p_{\beta}, \beta, q) = J(\beta). $$

We differentiate with respect to $\beta$ both sides of this identity:

$$ D_{\beta} J = D_{p} L \circ D_{\beta} p_{\beta} + D_{\beta} L. \quad (24) $$
The approach consists in choosing an adjoint variable $q$ such that $D_pL = 0$, which circumvents the difficulty to explicitly estimate $D_βp_β$. Let $\tilde{p}$ denote a variation in the variable $p$. It holds that

$$ 
\begin{aligned}
D_p L \tilde{p} &= - \int_{\Gamma_{out}} \frac{\partial \tilde{p}}{\partial n} + \int_{\Omega} \nabla \tilde{p} \cdot \nabla q - \int_{\Gamma_{out}} \frac{\partial \tilde{p}}{\partial n} q - \int_{\Gamma_{up}} \frac{\partial \tilde{p}}{\partial n} q \\
&= - \int_{\Gamma_{out}} \frac{\partial \tilde{p}}{\partial n} (q+1) + \int_{\Omega} (-\Delta q) \tilde{p} + \int_{\Gamma} \frac{\partial q}{\partial n} \tilde{p} - \int_{\Gamma_{up}} \frac{\partial \tilde{p}}{\partial n} q.
\end{aligned}
$$

Since $\tilde{p}$ is a variation of $p$, which identically vanishes on $\Gamma_0$, the boundary integral over $\Gamma = \partial\Omega$ reduces to an integral over $\Gamma_w$. The adjoint problem is designed in order to make the previous expression vanish:

$$ -\Delta q = 0 \text{ in } \Omega, q = -1 \text{ on } \Gamma_{out}, \frac{\partial q}{\partial n} = 0 \text{ on } \Gamma_w, q = 0 \text{ on } \Gamma_{up}, $$

which is exactly (22), up to a change in sign in the Dirichlet boundary condition. Let $q$ be the solution to this adjoint problem, so that $D_pL(p_β, β, q) = 0$. From (29), it comes

$$ D_β J = 0 + D_β L(p_β, β, q), $$

with

$$ D_β J \tilde{\beta} = \int_{\Omega} q \nabla \cdot \tilde{\beta} = - \int_{\Omega} \tilde{\beta} \cdot \nabla q + \int_{\Gamma} q \tilde{\beta} \cdot n. $$

Since $\tilde{\beta} = 0$ on $\Gamma_{out}$, $q = 0$ on $\Gamma_{up}$, and $U \cdot n = 0$ on $\Gamma_w$, the boundary term vanishes, which yields

$$ \nabla \beta J = -\nabla q, $$

which ends the proof. $\square$

<figure><img src="image_25.png" /></figure>

Figure 13 represents $-\nabla q$. The fact that it is the gradient of $J(β)$, that is the outflow through the exit, means that, at any point $x$ of the saturated domain, if a person located around $x$ change their desired velocity by a vector $β$, the effect upon the flux will be proportional to $-\nabla q \cdot \beta$. As shown in the figure, the gradient of $J(β)$ is strikingly similar to a desired velocity field, which means that any local variation of the desired velocity in its own direction will lead to an increase of the flux through the door, which corresponds to a *Faster is Faster* effect.

Let us now show that the situation is fully different in the microscopic setting. As previously, we shall investigate the dependence of some observable quantifying the evacuation upon variations of individual desired velocities of people upstream the exit. We aim here at determining whether it might be possible that some agents would improve the evacuation efficiency by reducing their pace. In the discrete situation, the instantaneous flux through the exit is not a relevant marker of egress efficiency. We shall rather consider the velocity of a particular agent close to the exit. We consider a highly congested situation like the one represented in Fig. 12. We propose to formulate the FiS issue in the following terms: consider a person (i.e. a disk) of index $i$ on the front of the cluster, i.e. close to the exit (see Fig. 12). We denote by $n_i$ the direction of their desired velocity. Now consider that any other individual $j$ has the ability to slightly change their speed, i.e. by changing their desired velocity to $U_j + \beta_j$. The FiS issue amounts to investigate whether taking $\beta_j$ in the direction of $U_j$ for some individuals might *decrease* the velocity of $i$ in the desired direction, thereby decreasing the flow near the exit. To formalize this question, we first consider the actual velocity field associated to the desired velocity field $U$ at some instant, i.e. the couple $(u, p)$ solves System (16). As in Proposition 25, we eliminate all those rows of $B$ that do not correspond to active constraints. We keep the same notation for the reduced $B$ and the associated reduced vector $p$ of Lagrange multipliers. The system writes

$$ 
\left| 
\begin{array}{l}
u + B^* p = U, \\
Bu = 0,
\end{array}
\right. 
\qquad (25) $$
<figure><img src="image_26.png" /><figcaption>**Figure 12.** Definition of $i$ and $n_i$.</figcaption></figure>

where all pressures are positive. The velocity can be eliminated, which yields

$$BB^* p = BU,$$

which is the discrete counterpart of the macroscopic Poisson problem ($-\Delta p = -\nabla \cdot U$).

We now consider the velocity as a function $u_\beta$ of $\beta \in \mathbb{R}^{2N}$, i.e. the desired velocity of individual $j$ is $U_j + \beta_j$. Admitting that the network of active contacts remains the same for small variations of $\beta$, it holds that

$$BB^* p_\beta = BU + B\beta, \quad u_\beta = U + \beta - B^* p_\beta.$$

The objective function we are interested in is the actual velocity of $i$ in its desired direction:

$$u_{\beta} \cdot n_{i} = U_{i} \cdot n_{i} + \beta \cdot n_{i} - B^{*} p_{\beta} \cdot n_{i}, \quad n_{i} = (0, \dots, 0, U_{i} / |U_{i}|, 0, \dots, 0). \qquad (26)$$

Since person $i$ is the closest to the exit, with no one in front, any change of $i$'s desired velocity shall affect their actual velocity in the same direction. We shall therefore disregard variations of speed for $i$, so that the first two terms above are constant, and we define the objective function as

$$J(\beta) = -B^* p_\beta \cdot n_i. \qquad (27)$$

**Proposition 31.** Let $\beta \mapsto J(\beta)$ be the flux functional defined by (27). The gradient of $J$ with respect to $\beta$ is $-B^* q$, where $q$ is the solution to the adjoint problem

$$BB^* q = B n_i. \qquad (28)$$

**Proof.** Like in the macroscopic setting, the core of the approach relies on a so-called Lagrangian, that is a function of the state variable $p$, the control variable $\beta$, and a dual variable $q$. It is defined as the sum of the objective function (velocity of $i$ along the desired direction) expressed in the state variable $p$ (uncoupled from the control variable $\beta$), and a weak expression of the state equation $BB^* p_\beta = BU + B\beta$:

$$L(p, \beta, q) = -B^* p \cdot n_i + (BB^* p - BU - B\beta) \cdot q.$$

As a consequence, for any $\beta$ and any $q$,

$$L(p_\beta, \beta, q) = J(\beta).$$

We differentiate with respect to $\beta$ both sides of this identity:

$$D_{\beta} J = D_{p} L \circ D_{\beta} p_{\beta} + D_{\beta} L. \qquad (29)$$
The approach consists again in choosing an adjoint variable $q$ such that $D_p L = 0$, which circumvents the difficulty to explicitly estimate $D_\beta p_\beta$. Let $\tilde{p}$ denote a variation in the variable $p$. It holds that

$$D_p L \tilde{p} = -B^* \tilde{p} \cdot n_i + (BB^* \tilde{p}) \cdot q = (-B n_i + BB^* q) \cdot \tilde{p}.$$

The adjoint problem is designed in order to make the previous expression vanish, it therefore reads

$$BB^* q = B n_i. \tag{30}$$

On the other hand, it holds that

$$D_\beta L \tilde{\beta} = -B \tilde{\beta} \cdot q = -B \star q \cdot \tilde{\beta}.$$

We finally obtain

$$\nabla J(\tilde{\beta}) = -B^* q$$

where $q$ is the solution to the adjoint problem (30). □

This expression of the gradient makes it possible to investigate the Faster is Slower effect. Figure 13 represents a congested situation, the black disk corresponds to $i$, the agent upfront. The black arrows represent to $\nabla J$, i.e. $-B^* q$, where $q$ is the solution to the adjoint problem. This field is the straight discrete counterpart to the continuous field represented on the top of the figure, yet it presents very different characteristics. Contrary to the macroscopic field, some arrows indeed point upstream the exit. For all such agents (represented in blue in the figure), changing their velocity against their desired direction (i.e. slowing down) will increase the speed of $i$, i.e. it will speed up the evacuation process.

## 5. Abstract considerations

We propose here an abstract setting to handle the macroscopic and microscopic saddle-point formulations (systems (12) and (16), respectively), and we give an abstract criterium to distinguish those situations which lead to a maximum principle from those which do not.

### Abstract saddle-point formulations

The abstract setting is represented in Fig. 14 (middle). The primal space $V$ is $L^2(\emptyset)$ in the macroscopic setting, and $\mathbb{R}^{2N}$ in the microscopic one. The figure represents the different cones which can be introduced. The left-hand side of the figure represents $C$, the feasible cone, and its polar

$$C^o = \{v \in V, \langle v | u \rangle \le 0 \quad \forall u \in C\},$$

in $(\ker B)^\perp$, so that $\ker B$ is reduced to single point, which is $C \cap C^o$, in this representation.

**Proposition 32.** Let $V$ and $\Lambda$ be Hilbert spaces, $B \in \mathcal{L}(V, \Lambda)$, $\Lambda_+$ a closed convex cone in $\Lambda$. We define its polar cone as

$$\Lambda_+^o = \{q \in \Lambda, \langle q | p \rangle \le 0 \quad \forall p \in \Lambda_+\}.$$

For a given $U \in V$, we consider the problem which consists in minimizing

$$v \mapsto J(v) = \frac{1}{2} |v - U|^2,$$

over $C \subset V$ defined by

$$C = \{v \in V, \langle Bv | q \rangle \le 0 \quad \forall q \in \Lambda_+\} = B^{-1}(\Lambda_+^o).$$

There exists a unique minimizer for $J$ over $C$ and, in any of the following situations:

(i) $\Lambda_+$ is spanned by a finite number of vectors, i.e. $\Lambda_+ = \text{co}(g_1, \dots, g_p)$ (conic convex hull).

(ii) $B$ is surjective.
<figure><img src="image_27.png" /><figcaption>**Figure 13.** Solution to the adjoint problem, in the macroscopic (top) and microscopic (bottom) settings.</figcaption></figure>

the problem admits a variational formulation in the following sense: $u$ is the constrained minimizer if and only if there exists $p$ such that

$$ \left| \begin{array}{ll} u + B^* p = U & \\ Bu \in \Lambda_+^\circ & \\ p \in \Lambda_+ & \\ \langle p | Bu \rangle = 0. & \end{array} \right. \qquad (31) $$

**Proof.** The Hilbert space $V$ is decomposed into the sum of $C$ and its polar cone $C^\circ$, i.e. any $U \in V$ decomposes uniquely as the sum

$$ U = u + u^{\circ}, \ u \in C, \ u^{\circ} \in C^{\circ}, \ \langle u | u^{\circ} \rangle = 0, $$

which can be seen as a unilateral version of the direct sum of two orthogonal subspaces (see [28]). Let us prove now that, in both cases (i) and (ii), $u^o$ can be written $B^*p$ with $p \in \Lambda_+$.

Case (i): If $\Lambda_+$ is the closed convex conic hull of a finite number of vectors, i.e. if

$$ \Lambda_+ = \text{co}(g_1, \dots, g_p) = \{ \sum_i \lambda_i g_i, \text{ with } \lambda_i \ge 0 \quad \forall i = 1, \dots, p \}, $$
<figure><img src="image_28.png" /><figcaption>**Figure 14.** Abstract setting.</figcaption></figure>

then C can be written

$$C = \{v \in V, \langle v | B^* g_i \rangle \le 0 \quad \forall i = 1, \dots, p\},$$

and the identity

$$C^\circ = \co(B^* g_1, \dots, B^* g_p) = \{ \sum_i \lambda_i B^* g_i, \text{ with } \lambda_i \ge 0 \quad \forall i = 1, \dots, p \}$$

is a straight consequence of Farkas' Lemma. Note that this lemma is usually expressed in a finite dimensional setting, but it can be easily checked that it generalizes to infinite dimensional spaces as soon as the number of unilateral constraints is finite (the cone has a finite number of facets). The previous identity implies that $u^\circ \in C^\circ$ can be written $u^\circ = B^*p$, with $p \in \lambda_+$, which settles case (i).

Case (ii): In all generality it holds that

$$C^\circ = \overline{B^* \Lambda_+}$$

Indeed, $C^\circ = \overline{B^* \Lambda_+}$ is trivial, and if the inclusion is strict, there exists $v \in \overline{B^* \Lambda_+}$, with $v \notin C^\circ$. By Hahn-Banach Theorem, there exist $w \in V$ and $\alpha \in \mathbb{R}$ such that

$$\langle B^* q | w \rangle \le \alpha < \langle v | w \rangle \quad \forall q \in \Lambda_+.$$

Since $\Lambda_+$ is a cone, the first inequality implies that $\langle B^* q | w \rangle \le 0$ for any $q \in \Lambda_+$, so that $w \in C$, and $\alpha \ge 0$. We then have $\langle v | w \rangle > 0$, with $v \in C^\circ$ and $w \in C$, which is impossible.

As a consequence, if $B^* \Lambda_+$ is closed, then $C^\circ = B^* \Lambda_+$. Let us prove that $B^* \Lambda_+$ is indeed closed. If $B$ is surjective between the Hilbert spaces $V$ and $\Lambda$, then there exists $\gamma > 0$ such that (continuous inf-sup condition)

$$\inf_{\Lambda \subset V} \sup_{V} \frac{\langle Bv | q \rangle}{|v| |q|} \ge \gamma \Rightarrow |B^* q| \ge \gamma |q| \quad \forall q \in \Lambda.$$
As a consequence, if a sequence $B^*q_n$ converges in $V$, with $q_n \in \Lambda_+$, then $q_n$ is a Cauchy sequence that converges to $q \in \Lambda_+$, and the previous limit is $B^*q \in B^*\Lambda_+$, which proves that $B^*\Lambda_+$ is closed, and thereby ends the proof. □

**Remark 33.** In the case of equality affine constraints, assuming that $B(V)$ is closed is sufficient to ensure existence of a Lagrange multiplier, whereas surjectivity ensures its uniqueness. In the present situation, with inequality constraints, assuming surjectivity, as we did, is essential to ensure existence. Indeed, an operator with closed range maps a closed convex cone to a convex cone which is not necessarily closed: if $B^*(\Lambda)$ is closed (which is equivalent to say that $B$ has a closed range), it does not imply that $B^*(\Lambda_+)$ is closed, even in the finite dimensional setting, we refer to [21] for a counter-example.

*Abstract maximal principle*

The difference between the macroscopic and the microscopic models can be described within the abstract setting presented above. In the macroscopic case, if one considers a desired velocity that is concentrating, i.e. $-\nabla \cdot U \ge 0$, it writes $BU \in -C$ (the operator $B$ is $-\nabla \cdot$). The fact that $-\Delta p = -\nabla \cdot U \ge 0$ (which writes $BB^*p = BU$), admits a solution that is non-negative means that there exists $p \in \Lambda_+$ such that $U - B^*p$ is in $kerB$. In other words, the maximal principle takes the abstract form

$$-C \subset C^o + \ker B,$$

which is equivalent to $-C \subset B^*(\Lambda_+) + \ker B$. Fig. 14 (middle), which gives a representation in $(\ker B)^\perp$, correspond to this situation : $-C$ appears as included in $C^o$, and the vector $U'$ (in green) corresponds to a concentrating field (i.e. with nonpositive divergence).

The microscopic situation corresponds to the situation represented in the bottom. The inclusion above no longer holds :

$$-C \not\subset B^*(\Lambda_+) + \ker B.$$

There exists vectors $U'$ in $-C$, i.e. concentrating fields like the one represented in Fig. 7 (top-right), that are not in $C^o + \ker B$ so that there exists no pressure fields in $\Lambda_+ = \mathbb{R}_+^{N_c}$ such that $BB^*p = BU$, which expresses the pathology of this Laplacian $BB^*$.

## 6. Conclusion

We have recalled in this paper the tight links between continuous and discrete Laplacians, in contexts where the discrete Laplacian combines a phenomenological law of the Ohm's type together with a conservation principle. In the context of crowd motions (or granular flows), accounting for the non-overlapping constraint between solid grains can be done by duality, which leads to a unilateral discrete Darcy problem. Like in the continuous setting, eliminating the velocity leads to a Poisson-like problem for a pressure field defined on the dual graph of the contact network. The discrete Laplacian $BB^*$ is the straight counterpart of $-\nabla \cdot \nabla = -\Delta$ in the one-dimensional situation (people / grains in a row), but the situation is fully different in higher dimensions. In particular, the matrix $BB^*$ does not verify any maximum principle in general, and it does vanish on constant fields. We have investigated how these very defects make it possible to recover effects that are observed in reality, whereas in the continuous setting, with a standard Laplacian verifying the maximum principle, these effects are not reproduced. These considerations also rule out the possibility to properly obtain the continuous model as a many-body limit of the discrete one, as illustrated for example by Figs. 11 and 13, which enlight the deep discrepancies between both scales, although the corresponding models are built upon the very same principles.
## Conflicts of interest

The authors declare no competing financial interest.

## Acknowledgments

The author would like to thank S. Faure for his contributions in numerical computations, and R. Glowinski for stimulating discussions over the last decades, which lead to these developments.

## References

[1] B. Maury, *The Respiratory System in Equations*, MS&A, Springer Milano, 2013.

[2] C. A. Whitfield, P. Latimer, A. Horsley, J. M. Wild, G. J. Collier, O. E. Jensen, "Spectral graph theory efficiently characterizes ventilation heterogeneity in lung airway networks", *J. R. Soc. Interface* **17** (2020).

[3] J. Droniou, "Finite volume schemes for diffusion equations: Introduction to and review of modern methods", *Mathematical Models and Methods in Applied Sciences* **24** (2014), no. 8, p. 1453-1455.

[4] J. Ferreira, M. Grassi, E. Gudiño, P. de Oliveira, "A new look to non-Fickian diffusion", *Applied Mathematical Modelling* **39** (2015), no. 1, p. 194-204.

[5] J.-L. Lions, E. Magenes, *Non-Homogeneous Boundary Value Problems and Applications*, 1 ed., vol. 1, Grundlehren der mathematischen Wissenschaften, no. 181, Springer Berlin, Heidelberg, 1972.

[6] E. Rapetti, A. Bossavit, "Whitney Forms of Higher Degree", *SIAM Journal on Numerical Analysis* **47** (2009), no. 3, p. 2369-2386, https://doi.org/10.1137/070705489, https://arxiv.org/abs/https://doi.org/10.1137/070705489.

[7] J. L. Doob, *Classical Potential Theory and Its Probabilistic Counterpart*, Grundlehren der mathematischen Wissenschaften, Springer New York, NY, 1984.

[8] D. Gilbarg, N. S. Trudinger, *Elliptic Partial Differential Equations of Second Order*, 2 ed., Grundlehren der mathematischen Wissenschaften, vol. 224, Springer Berlin, Heidelberg, 2001.

[9] I. Karatzas, S. Shreve, *Brownian Motion and Stochastic Calculus*, Graduate Texts in Mathematics, vol. 113, Springer, 1991.

[10] B. Gustafsson, "Lectures on Balayage", Report series, no. 7, p. 17-63, Report series **7**, Sirkka-Liisa Eriksson, 2002.

[11] M. E. Muller, "Some Continuous Monte Carlo Methods for the Dirichlet Problem", *The Annals of Mathematical Statistics* **27** (1956), no. 3, p. 569-589, https://doi.org/10.1214/aoms/1177728169.

[12] R. Lyons, Y. Peres, *Probability on Trees and Networks*, Cambridge University Press, 2017.

[13] R. Glowinski, B. Maury, "Fluid-Particle Flow: a Symmetric Formulation", *C. R. Acad. Sci. Paris, Série I* **324** (1997), p. 1079-1084.

[14] R. Glowinski, T. Pan, T. Hesla, D. Joseph, "A distributed Lagrange multiplier/fictitious domain method for particulate flows", *International Journal of Multiphase Flow* **25** (1999), no. 5, p. 755-794.

[15] B. Maury, "A time-stepping scheme for inelastic collisions", *Numerische Mathematik* **102** (2006), p. 649-679.

[16] B. Maury, A. Preux, "Pressureless Euler equations with maximal density constraint : a time-splitting scheme", *Topological Optimization and Optimal Transport: In the Applied Sciences*, chap. 13, Topological Optimization and Optimal Transport: In the Applied Sciences, Berlin, Boston: De Gruyter, 2017.

[17] P. Degond, P. Minakowski, L. Navoret, E. Zatorska, "Finite volume approximations of the Euler system with variable congestion", *Computers & Fluids* **169** (2018), p. 23-39, https://www.sciencedirect.com/science/article/pii/S0045793017303341.

[18] P. C. Bianchini, R., "Soft congestion approximation to the one-dimensional constrained Euler equations", *Nonlinearity* **34** (2021), no. 10.

[19] B. Maury, A. Roudneff-Chupin, F. Santambrogio, "A macroscopic Crowd Motion Model of the gradient-flow type", *Mathematical Models and Methods in Applied Sciences* **20** (2010), no. 10, p. 1787-1821.

[20] R. Trémolières, J.-L. Lions, R. Glowinski, *Numerical Analysis of Variational Inequalities*, Studies in Mathematics and its Applications, North Holland, 2000.

[21] F. Bourdin, B. Maury, "Multibody and macroscopic impact laws: a Convex Analysis standpoint", *Trails in Kinetic Theory: Foundational Aspects and Numerical Methods*, p. 97-139, *Trails in Kinetic Theory: Foundational Aspects and Numerical Methods*, Springer International Publishing, Cham, 2021.

[22] B. Maury, J. Venel, "A discrete contact model for crowd motion", *ESAIM: M2AN* **45** (2011), no. 1, p. 145-168.

[23] E. Cepolina, N. Tyler, "Understanding Capacity Drop for designing pedestrian environments", in *Proceedings of Walk 21: Everyday Walking Culture*, 2005.

[24] C. D. B. Andreianov, M. Rosini, "Crowd dynamics and conservation laws with nonlocal constraints and capacity drop", *Math. Models Methods Appl. Sci.* (2014), p. 2685-2722.
[25] B. Andreianov, C. Donadello, U. Razafison, M. D. Rosini, "Qualitative behaviour and numerical approximation of solutions to conservation laws with non-local point constraints on the flux and modeling of crowd dynamics at the bottlenecks", *M2AN* **50** (2016), no. 5, p. 1269-1287.

[26] C. Gershenson, D. Helbing, "When slower is faster", *Complexity* **21** (2015), p. 9-15.

[27] B. Maury, S. Faure, *Crowds in Equations*, London: World Scientific Publishing Europe Ltd, London, 2018.

[28] J.-J. Moreau, "Décomposition orthogonale d'un espace Hilbertien selon deux cônes mutuellement polaires", *C. R. Acad. Sci. Paris* **255** (1962), p. 238-240.

[29] A. Bossavit, "Mixed finite elements and the complex of Whitney forms", The mathematics of finite elements and applications, vol. VI, The mathematics of finite elements and applications, Academic Press, London, 1988.

[30] K. Yuan, V. L. Knoop, S. P. Hoogendoorn, "Capacity Drop: Relationship Between Speed in Congestion and the Queue Discharge Rate", *Transportation Research Record* **2491** (2015), no. 1, p. 72-80, https://doi.org/10.3141/2491-08, https://arxiv.org/abs/https://doi.org/10.3141/2491-08.

[31] S. Faure, B. Maury, "Crowd motion from the granular standpoint", *Mathematical Models and Methods in Applied Sciences* **25** (2015), no. 03, p. 463-493.

[32] B. Maury, A. Roudneff-Chupin, F. Santambrogio, J. Venel, "Handling Congestion in Crowd Motion Modeling", *Networks and Heterogeneous Media* **6** (2011), no. 3, p. 485-519.