# Nonlocal-Interaction Equation on Graphs: Gradient Flow Structure and Continuum Limit

Antonio Esposito, Francesco S. Patacchini, André Schlichting, Dejan Slepčev

► To cite this version:

Antonio Esposito, Francesco S. Patacchini, André Schlichting, Dejan Slepčev. Nonlocal-Interaction Equation on Graphs: Gradient Flow Structure and Continuum Limit. Archive for Rational Mechanics and Analysis, 2021, 240 (2), pp.699-760. 10.1007/s00205-021-01631-w. hal-03227655

HAL Id: hal-03227655

https://ifp.hal.science/hal-03227655

Submitted on 17 May 2021

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Nonlocal-Interaction Equation on Graphs:
## Gradient Flow Structure and Continuum Limit

ANTONIO ESPOSITO, FRANCESCO S. PATACCHINI, ANDRÉ SCHLICHTING & DEJAN SLEPCEV

Communicated by A. FIGALLI

## Abstract

We consider dynamics driven by interaction energies on graphs. We introduce graph analogues of the continuum nonlocal-interaction equation and interpret them as gradient flows with respect to a graph Wasserstein distance. The particular Wasserstein distance we consider arises from the graph analogue of the Benamou-Brenier formulation where the graph continuity equation uses an upwind interpolation to define the density along the edges. While this approach has both theoretical and computational advantages, the resulting distance is only a quasi-metric. We investigate this quasi-metric both on graphs and on more general structures where the set of “vertices” is an arbitrary positive measure. We call the resulting gradient flow of the nonlocal-interaction energy the nonlocal nonlocal-interaction equation (NL²IE). We develop the existence theory for the solutions of the NL²IE as curves of maximal slope with respect to the upwind Wasserstein quasi-metric. Furthermore, we show that the solutions of the NL²IE on graphs converge as the empirical measures of the set of vertices converge weakly, which establishes a valuable discrete-to-continuum convergence result.

## Contents

<table><tr><td>1. Introduction</td><td>701</td></tr><tr><td>1.1. Graph Setting with General Interactions</td><td>701</td></tr><tr><td>1.2. General Setting for Vertices in Euclidean Space</td><td>704</td></tr><tr><td>1.3. Relation to the Numerical Finite-Volume Upwind Scheme</td><td>708</td></tr><tr><td>1.4. Comparison with Other Discrete Metrics and Gradient Structures</td><td>709</td></tr><tr><td>1.5. Connections to Machine Learning</td><td>712</td></tr><tr><td>2. Nonlocal Continuity Equation and Upwind Transportation Metric</td><td>714</td></tr><tr><td>2.1. Weight Function</td><td>714</td></tr><tr><td>2.2. Action</td><td>715</td></tr><tr><td>2.3. Nonlocal Continuity Equation</td><td>720</td></tr><tr><td>2.4. Nonlocal Upwind Transportation Quasi-Metric</td><td>725</td></tr></table>
<table><tr><td>3. Nonlocal Nonlocal-Interaction Equation</td><td>732</td></tr><tr><td>3.1. Identification of the Gradient in Finsler Geometry</td><td>734</td></tr><tr><td>3.2. Variational Characterization for the Nonlocal Nonlocal-Interaction Equation</td><td>738</td></tr><tr><td>3.3. The Chain Rule and Proof of Theorem 3.9</td><td>742</td></tr><tr><td>3.4. Stability and Existence of Weak Solutions</td><td>747</td></tr><tr><td>3.5. Discussion of the Local Limit</td><td>753</td></tr><tr><td>Acknowledgement</td><td>754</td></tr><tr><td>Appendix</td><td>755</td></tr><tr><td>References</td><td>757</td></tr></table>

Notation

We list here some symbols used throughout the paper.

* $\mathcal{M}(A)$ is the set of Borel measures on $A \subseteq \mathbb{R}^d$.

* $\mathcal{M}^+(A)$ is the set of non-negative Borel measures on $A$.

* $\mathcal{P}(A) \subset \mathcal{M}^+(A)$ is the set of Borel probability measures on $A$.

* $\mathcal{P}_2(A) \subseteq \mathcal{P}(A)$ stands for the elements of $\mathcal{P}(A)$ with finite second moment, that is,

$$M_2(\rho) := \int_A |x|^2 \, d\rho(x) < \infty.$$

* $C_b(A)$ is the set of bounded continuous functions from $A$ to $\mathbb{R}$.

* $a_+ := \max\{0, a\}$ and $a_- := (-a)_+$ are the positive and negative parts of $a \in \mathbb{R}$.

* $\mu \in \mathcal{M}^+(\mathbb{R}^d)$ sets the underlying geometry of the state space; it is sometimes referred to as base measure.

* $\rho \in \mathcal{P}(\mathbb{R}^d)$ denotes a configuration; the natural setting is that $\text{supp} \rho \subseteq \text{supp} \mu$, although we allow for general supports as needed for stability results.

* $\eta: \{(x, y) \in \mathbb{R}^d \times \mathbb{R}^d : x \neq y\} \to [0, \infty)$ is the edge weight function.

* $G = \{(x, y) \in \mathbb{R}^d \times \mathbb{R}^d : x \neq y, \eta(x, y) > 0\}$ are the edges.

* $\rho_1 \otimes \rho_2 \in \mathcal{M}^+(G)$ is the product measure of $\rho_1, \rho_2 \in \mathcal{M}^+(\mathbb{R}^d)$ restricted to $G$.

* $\gamma_1 = \rho \otimes \mu$ and $\gamma_2 = \mu \otimes \rho$.

* $\mathcal{V}^{\text{as}}(G)$ is the set of antisymmetric graph vector fields on $G$, defined in (1.6).

* $\bar{\nabla}f$ is the nonlocal gradient of a function $f: \mathbb{R}^d \to \mathbb{R}$, while $\bar{\nabla} \cdot j$ is the nonlocal divergence of a measure-valued flux $j \in \mathcal{M}(G)$; see Definition 2.7.

* $A$ stands for the action functional; see Definition 2.3.

* $\mathcal{T}$ denotes the nonlocal transportation quasi-metric; see (2.22).

* $\text{CE}_T(\rho_0, \rho_1)$ denotes the set of paths (solutions to the nonlocal continuity equation for densities (1.7) or measures (2.12)) on the time interval $[0, T]$ connecting two measures $\rho_0, \rho_1 \in \mathcal{P}(\mathbb{R}^d)$; we set $\text{CE} := \text{CE}_1$.

Let us also specify the notions of *narrow convergence* and *convolution*. A sequence $(\rho^n)_n \subset \mathcal{M}(A)$ is said to converge narrowly to $\rho \in \mathcal{M}(A)$, in which case we write $\rho^n \to \rho$, provided that

$$\forall f \in C_b(A), \quad \int_A f \, d\rho^n \to \int_A f \, d\rho \quad \text{as } n \to \infty.$$
Given a function $f: A \times A \to \mathbb{R}$ and $\rho \in \mathcal{M}(A)$, we write $f * \rho$ the convolution of $f$ and $\rho$, that is,

$$f * \rho(x) = \int_A f(x, y) \, d\rho(y) \quad \text{for any } x \in A \text{ such that the right-hand side exists.}$$

# 1. Introduction

We investigate dynamics driven by interaction energies on graphs, and their continuum limits. We interpret the relevant dynamics as gradient flows of the interaction energy with respect to a particular graph analogue of the Wasserstein distance. We prove the convergence of the dynamics on finite graphs to a continuum dynamics as the number of vertices goes to infinity. To do this we create a unified setup where the continuum and the discrete dynamics are both seen as particular instances of the gradient flow of the same energy, with respect to a nonlocal Wasserstein quasi-metric whose state space is adapted to the configuration space considered.

Let us first introduce the problem on finite graphs where it is the simplest to describe.

## 1.1. Graph Setting with General Interactions

Consider an undirected graph with vertices $X = \{x_1, \dots, x_n\}$ and edge weights $w_{x,y} \ge 0$, satisfying $w_{x,y} = w_{y,x}$ for all $x, y \in X$. Although technically not necessary, we impose the natural requirement that $w_{x,x} = 0$. The interaction potential is a symmetric function $K: X \times X \to \mathbb{R}$, while the external potential is denoted $P: X \to \mathbb{R}$. We consider a “mass” distribution $\rho: X \to [0, \infty)$, and we require $\sum_{x \in X} \rho_x = 1$. The total energy $\mathcal{E}_X: \mathcal{P}(X) \to \mathbb{R}$ is a combination of the interaction energy $\mathcal{E}_I$ and the potential energy $\mathcal{E}_P$:

$$\mathcal{E}_X(\rho) = \mathcal{E}_I(\rho) + \mathcal{E}_P(\rho) = \frac{1}{2} \sum_{x \in X} \sum_{y \in X} K_{x,y} \rho_x \rho_y + \sum_{x \in X} P_x \rho_x. \quad (1.1)$$

The gradient descent of $\mathcal{E}_X$ that we study is described by the following system of ODE for the mass distribution:

$$\frac{d\rho_x}{dt} = -\frac{1}{2} \sum_{y \in X} (j_{x,y} - j_{y,x}) w_{x,y}, \quad (1.2)$$

$$j_{x,y} = \frac{1}{n} (\rho_x(v_{x,y})_+ - \rho_y(v_{x,y})_-), \quad (1.3)$$

$$v_{x,y} = - \sum_{z \in X} \rho_z (K_{y,z} - K_{x,z}) - (P_y - P_x). \quad (1.4)$$

The quantities $v: X \times X \to \mathbb{R}$ and $j: X \times X \to \mathbb{R}$ are defined on edges and model the graph analogues of velocity and flux. An evolution by such system is illustrated on Fig. 1. The system (1.2)–(1.4) is the gradient flow of the energy $\mathcal{E}_X$ with respect
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** A solution of the nonlocal-interaction equation on graphs driven by the energy (1.1). We consider a graph based on 240 sample points $X$ from a 2D two-moon data set. The connectivity distance is $\varepsilon = 0.7$. The edge weights are $w_{x,y} = \exp(-6|x-y|^2)$ if $|x-y| \le \varepsilon$ and zero otherwise. The interaction potential is $K_{x,y} = 1 - \exp(-d(x,y)^2/10)$, where $d(x,y)$ is the graph distance between vertices $x$ and $y$ of $X$ with edge weights $1/w_{x,y}$, and the external potential is $P = 0$. The solution, starting from a uniform distribution, is shown at time $t = 60$. Brighter color indicates more mass</figcaption></figure>

to a new graph equivalent of the Wasserstein metric. The concept of Wasserstein metrics on finite graphs were introduced independently by CHOW et al. [14], MAAS [36], and MIELKE [37,38]. All of the approaches rely on graph analogues of the continuity equation to describe the paths in the configuration space. On graphs the mass is distributed over the vertices and is exchanged over the edges. Hence, the analogues of the vector field and the flux are defined over the edges. However, the flux should be the product of the velocity (an edge-based quantity) by the density (a vertex-based quantity). Thus, one has to interpolate the densities at vertices to define the density (and hence the flux) along the edges. The choice of interpolation is not unique, and has important ramifications.

While the overall structure of our setup is derived from one in [36], which we recall in Section 1.4; the form of the interpolation used is related to the upwind interpolation used in [14] and is almost identical to one in [13]. While in [14] the authors considered only the direction of the flux due to the potential energy to determine which density to use on the edges, in our case the density chosen depends on the total velocity and we furthermore include the interaction term which itself depends on the configuration. In particular, we use an upwind interpolation based on the total velocity. In the context of graph Wasserstein distance, such interpolation was first used by CHEN et al. [13].

The “velocities” $v$ we consider can be assumed to be antisymmetric: $v_{x,y} = -v_{y,x}$ for all $x, y \in X$. In the graph setting, which we normalize in order to consider
limit $n \to \infty$, the continuity equation with upwind interpolation is obtained by combining (1.2) with the flux-velocity relation (1.3). Similarly to [36] and exactly as in [13], we define the graph Wasserstein distance by minimizing the action, which is the graph analogue of the kinetic energy:

$$A(\rho, v) = \frac{1}{n} \sum_{x \in X} \sum_{y \in X} (v_{x,y+})^2 w_{x,y} \rho_x.$$

As in [13,14,36,38], the graph Wasserstein distance is defined by adapting the Benamou-Brenier formula:

$$\mathcal{T}(\rho^0, \rho^1)^2 = \inf_{(\rho, v) \in \mathcal{CE}_X(\rho^0, \rho^1)} \int_0^1 A(\rho(t), v(t)) \, dt,$$

where $\mathcal{CE}_X(\rho^0, \rho^1)$ is the set of all paths (i.e., solutions of (1.2)-(1.3)) connecting $\rho^0$ and $\rho^1$.

It is important to observe that, in our setting, $\mathcal{T}$ is not symmetric (that is, $\mathcal{T}(\rho^0, \rho^1)$ is in general different from $\mathcal{T}(\rho^1, \rho^0)$). The reason for this is that in general, $A(\rho, v) \neq A(\rho, -v)$. Therefore the nonlocal Wasserstein distance which arises from the upwind interpolation is only a quasi-metric. The action $A(\rho, v)$ provides a Finsler structure to the tangent space, instead of the usual Riemannian structure. Formally the system (1.2)-(1.4) is the gradient flow of $\mathcal{E}_X$ with respect to this Finsler structure; we present a derivation of this fact in a more general setting in Section 3.1. The system is also the curve of steepest descent with respect to quasi-metric $\mathcal{T}$, which is the point of view we use to create rigorous theory in the general setting.

**Remark 1.1.** The well-posedness of (1.2)-(1.4) is a straightforward consequence of the Picard existence theorem. Namely, note that the simplex $1 \ge \rho_x \ge 0$, $\sum_{x \in X} \rho_x = 1$ is an invariant region of the dynamics and that on it the vector field (1.4) is Lipschitz continuous in $\rho_x$, $x \in X$.

**Remark 1.2.** One could consider other interpolations instead of the upwind one. In particular, if we considered an interpolation of the form $I(\rho_x, \rho_y)$ instead of the upwind one, the only change in the gradient flow would be that the velocity-flux relation (1.3) would become $j_{x,y} = \frac{1}{n} I(\rho_x, \rho_y) v_{x,y}$. We note that this can have major implications on the resulting dynamics. In particular, for the logarithmic interpolation, $I(r, s) = (r-s)/(\ln r - \ln s)$, or the geometric interpolation, $I(r, s) = \sqrt{rs}$, the resulting dynamics would never expand the support of the solutions, so even for repulsive potentials the mass may not spread throughout the domain. On the other hand, using the arithmetic interpolation, $I(r, s) = (r+s)/2$, would not work directly since the solutions may become negative. In this case additional technical steps, like a Lagrange multiplier as in [39], are necessary to obtain the evolution of a non-negative probability density. We use the more physical inspired upwind flux, which automatically ensures the positivity of the density.

Before we turn to the general setting we point out that the system (1.2)-(1.4) offers a new model of graph-based clustering, which is briefly discussed in Section 1.5.
## 1.2. General Setting for Vertices in Euclidean Space

Here we introduce the general framework for studies of interaction equations on families of graphs and their limits as the number of vertices $n$ goes to $\infty$. In particular, in the applications to machine learning which we briefly discuss in Section 1.5, the graphs considered are random samples of some underlying measure in Euclidean space, and the edge weights, as well as the interaction energy, depend on the positions of the vertices. The vertices are points in $\mathbb{R}^d$. The edges are given in terms of a non-negative symmetric weight function $\eta: \{(x, y) \in \mathbb{R}^d \times \mathbb{R}^d : x \neq y\} \to [0, \infty)$, which defines the set of edges as $G = \{(x, y) \in \mathbb{R}^d \times \mathbb{R}^d : x \neq y, \eta(x, y) > 0\}$. From the discrete setting, the set of vertices is replaced by the more general notion of a measure on $\mathbb{R}^d$; the discrete graphs with vertices $X = \{x_1, \dots, x_n\} \subset \mathbb{R}^d$ correspond to $\mu$ being the empirical measure of the set of points, $\mu = \frac{1}{n} \sum_{i=1}^n \delta_{x_i}$. The distribution of mass over the vertices is described by the measure $\rho \in \mathcal{P}(\mathbb{R}^d)$ and in most applications we consider $\text{supp} \rho \subseteq \text{supp} \mu$. However, in order to prove general stability results (e.g., Theorem 3.14), we need to allow that initially part of the support of $\rho$ is outside of the support of $\mu$; we think of such mass as outside of the domain specified by $\mu$. The mass starting outside of the support of $\mu$ can only flow into the support of $\mu$. Here we present the evolution assuming $\rho \ll \mu$, while in Sections 2 and 3 we present the setup in full generality. Furthermore, we denote by $\rho$ both the measure and its density with respect to $\mu$.

The evolution of interest is the gradient descent of the energy $\mathcal{E}: \mathcal{P}(\mathbb{R}^d) \to \mathbb{R}$ given by

$$ \mathcal{E}(\rho) = \frac{1}{2} \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} K(x, y) \, d\rho(x) \, d\rho(y) + \int_{\mathbb{R}^d} P(x) \, d\rho(x), \quad (1.5) $$

where $K: \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}$ is symmetric and $P: \mathbb{R}^d \to \mathbb{R}$. This energy generalizes (1.1) in terms of the configurations $\rho$ and specializes it in terms of the type of potentials $K$ and $P$ considered. In fact, from now on we omit the subscripts $X$ referring to the vertices (e.g. in the energy) since our general setting allows for distribution of mass outside of the support of $\mu$. The gradient flow we consider takes the form

$$ \begin{align} \partial_t \rho_t(x) &= - \int_{\mathbb{R}^d} j_t(x, y) \eta(x, y) \, d\mu(y) =: -(\bar{\nabla} \cdot j_t)(x), \tag{NL$^2$IE} \\ j_t(x, y) &= \rho_t(x) v_t(x, y)_+ - \rho_t(y) v_t(x, y)_-, \nonumber \\ v_t(x, y) &= -(K * \rho_t(y) - K * \rho_t(x) + P(y) - P(x)). \nonumber \end{align} $$

The system (NL$^2$IE) consists first of a nonlocal continuity equation, where the divergence $\bar{\nabla}$. is encoded with the graph structure described through $\mu$ and $\eta$ (see Definition 2.7). Secondly, it involves a mapping from velocity to flux, which in our case is the upwind flux and encodes the geometry of the gradient structure. Finally, the third equation identifies the driving velocity as the nonlocal gradient of the variation of the energy (1.5). Overall, we obtain that (NL$^2$IE) is the gradient flow of the energy $\mathcal{E}$ with respect to a generalization of the graph Wasserstein metric we now introduce.
**Nonlocal Continuity Equation** Let us set

$$ \mathcal{V}^{\text{as}}(G) := \{v: G \to \mathbb{R}: v(x, y) = -v(y, x) \text{ for all } (x, y) \in G\} \quad (1.6) $$

and call its elements *nonlocal (antisymmetric) vector fields* on $G$; for any pair $(x, y) \in G$ the value $v(x, y)$ can be regarded as a jump rate from $x$ to $y$. Let us fix a final time $T > 0$ throughout the paper and let a family $\{v_t\}_{t \in [0, T]} \subset \mathcal{V}^{\text{as}}(G)$ be given. In the case $\rho_t \ll \mu$ for all $t \in [0, T]$, it is possible to combine the first two equations in (NL$^2$IE) in order to arrive at the nonlocal continuity equation

$$ \partial_t \rho_t(x) + \int_{\mathbb{R}^d} (\rho_t(x)v_t(x, y)_+ - \rho_t(y)v_t(x, y)_-) \eta(x, y) \, d\mu(y) = 0, \mu\text{-a.e. } x \in \mathbb{R}^d. \quad (1.7) $$

For general curves $\rho: [0, T] \to \mathcal{P}(\mathbb{R}^d)$, it is necessary to consider the weak form of (1.7), which is discussed in Section 2.3.

We remark that the general setup we develop allows for the solution $\rho$ to develop atoms and persist even after the atoms have formed. Heuristic arguments and numerical experiments indicate that there are equations covered by our theory for which this is the case. For example, if $\mu$ is the Lebesgue measure on $\mathbb{R}$, $\rho_0$ the restriction of the Lebesgue measure to $[-0.5, 0.5]$, $K(x, y) = |x - y|$ and $\eta(x, y) = 1/(x - y)^2$, then the solutions develop delta mass concentrations at 0 in finite time. Understanding for which $K$ and $\eta$ solutions do develop finite time singularities is an interesting open problem.

We note that when defining the flux in (1.7) we define the density along edges to be the density at the source; analogously to an upwind numerical scheme. While, as we show, this leads to a convenient framework to consider the dynamics, it creates the difficulty that the resulting distance, that we are about to define, is not symmetric and is thus only a quasi-metric.

**Upwind Nonlocal Transportation Metric** We use the nonlocal continuity equation (1.7) to define a nonlocal Wasserstein quasi-distance in analogy to the Benamou-Brenier formulation [6] for the classical Kantorovich-Wasserstein distances [50]. That is, for two probability measures $\rho_0, \rho_1 \in \mathcal{P}_2(\mathbb{R}^d)$, let

$$ T_{\mu}(\rho_0, \rho_1)^2 := \inf_{(\rho, v) \in \text{CE}(\rho_0, \rho_1)} \left\{ \iint_G \int_0^1 \int_0^1 |v_t(x, y)_+|^2 \eta(x, y) \, d\rho_t(x) \, d\mu(y) \, dt \right\}, \quad (1.8) $$

where $\text{CE}(\rho_0, \rho_1)$ is the set of weak solutions $\rho$ to the nonlocal continuity equation (see Definition 2.14) on $[0, 1]$ with $\rho(0) = \rho_0$ and $\rho(1) = \rho_1$. We note that the notion of the nonlocal Wasserstein distance for measures on $\mathbb{R}^d$ was introduced by ERBAR [23], who used it to study the fractional heat equation. One difference is that the interpolation we consider is beyond the scope of [23]. Very recently [43] has extended the gradient flow viewpoint of the jump processes to generalized gradient structures driven by a broad class of internal energies.

Another difference is that here the measure $\mu$ plays an important role in how the action is measured and allows one to incorporate seamlessly both the continuum
case (e.g., $\mu$ is the Lebesgue measure on $\mathbb{R}^d$) and the graph case ($\mu$ is the empirical measure of the set of vertices).

The notions above are rigorously developed in Section 2, where we list the precise assumption (W) on the edge weight function $\eta$ and the joint assumptions (A1) and (A2) on $\eta$ and the underlying measure $\mu$. We then rigorously introduce the action (Definition 2.3), which is a nonlocal analogue of kinetic energy; we show its fundamental properties, in particular joint convexity (Lemma 2.12) and lower semicontinuity with respect to narrow convergence (Lemma 2.9). In Section 2.3 we rigorously introduce the nonlocal continuity equation in measure-valued flux form (2.12); we introduce the notion on all of $\mathbb{R}^d$ where $\mu$ does not initially play a role. The measure $\mu$ enters the framework by considering paths of finite action. Proposition 2.17 establishes an important compactness property of sequences of solutions. In Section 2.4 we turn our attention to the nonlocal Wasserstein quasi-metric based on the upwind interpolation, which we introduce in Definition 2.18. The compactness of solutions of the nonlocal continuity equation and the lower semicontinuity of the action imply the existence of (directed) geodesics (Proposition 2.20). We do not characterize the geodesics. Nevertheless we note that this is a interesting problem. A possible approach in this direction is via duality using nonlocal analogues of the Hamilton-Jacobi equations, similarly to how this problem was recently treated in the discrete setting in [28, 30]. Following the work of ERBAR [23] we show that the nonlocal Wasserstein quasi-metric generates a topology on the set of probability measures which is stronger than the $W_1$ topology (i.e., the Monge distance or the 1-Wasserstein metric). Analogously to [2] we show the equivalence between the paths of finite length with respect to the quasi-metric and the solutions of the nonlocal continuity equation with finite action (Proposition 2.20). The set of probability measures endowed with the quasi-metric $\mathcal{T}$ has a formal structure of a Finsler manifold, and parts of this structure can be described; in particular, in (2.27) we describe the tangent space at a given measure $\rho$ using the fluxes. We note that using fluxes, instead of velocities, is necessary since, because of the upwinding, the relation between the velocities and the tangent vectors is not linear (Proposition 2.26) and in particular not symmetric. For this reason the resulting gradient structure is also different to the large class of nonlinear, however still symmetric, flux-velocity relations considered in [43]. We conclude Section 2 by showing that, given a measure $\mu$, the finiteness of the action ensures that any path starting within the support of $\mu$ will remain within the support of $\mu$ (Proposition 2.28).

**Nonlocal Nonlocal-Interaction Equation** In Section 3 we develop the existence theory of the equation (NL$^2$IE) based on the interpretation as the gradient flow of $\mathcal{E}$ with respect to the quasi-metric $\mathcal{T}$ defined in (1.8). We begin by listing the precise conditions (K1)–(K3) on the interaction kernel $K$. We note that these are less restrictive than the typical conditions for the well-posedness of the standard nonlocal-interaction equation in Euclidean setting [2, 10].

Before we turn to the rigorous theory of weak solutions as curves of maximal slope on quasi-metric space, we discuss the gradient flow structure in a more geometric setting, namely the Finsler structure related to $\mathcal{T}$. Indeed, the action [formally given by the time integrand in (1.8), and rigorously defined by (2.4)] defines
a positively homogeneous norm (namely a Minkowski norm) on the tangent space. The Hessian of the square of the norm endows the tangent space at each measure with the formal structure of a Riemann manifold. We compute this Riemann metric in “Appendix A” under an absolute-continuity assumption. With this assumption, we show that (NL²IE) is the gradient flow of $\mathcal{E}$ with respect to the Finsler structure in Section 3.1. For simplicity, we consider $P \equiv 0$, since the extension to $P \neq 0$ is straightforward, as it is explained in Remark 3.2.

In Section 3.2 we develop the rigorous gradient descent formulation based on curves of maximal slope in the space of probability measures endowed with the quasi-metric $\mathcal{T}$. The theory of gradient flows in the spaces of probability measures endowed with the standard Wasserstein metric was developed in [2]. Here we extend it to the setting of quasi-metric spaces, endowed with the nonlocal Wasserstein distance. This requires several delicate arguments. We start by introducing the notions of one-sided strong upper gradient (Definition 3.12) and curves of maximal slope (Definition 3.8). We define the local slope $\mathcal{D}$ in (3.19) by using a heuristically derived gradient of the energy $\mathcal{E}$, and show, using a chain rule established in Proposition 3.10, that $\sqrt{\mathcal{D}}$ is a one-sided strong upper gradient for $\mathcal{E}$ with respect to $\mathcal{T}$. One of our main results is Theorem 3.9, which establishes the equivalence between curves of maximal slope and weak solutions of (NL²IE). In Section 3.4 we prove several important results. Namely Theorem 3.14 establishes that the De Giorgi functional $\mathcal{G}_{\mathcal{T}}$ is stable under variations of the base measure $\mu$ and of the solutions. A consequence of this result is the convergence of solutions of (NL²IE) on graphs defined on random samples of a measure to solutions of (NL²IE) corresponding to the full underlying measure (Remark 3.17). The proof of Theorem 3.14 relies on the lower semicontinuity of the local slope (Lemma 3.12) and the lower semicontinuity of the De Giorgi functional (3.13). Another important consequence is the existence of weak solutions of (NL²IE), which is proved in Theorem 3.15.

**Remark 1.3.** (Asymptotics) Describing the steady states and determining the long-time asymptotics of (NL²IE) are natural and important problems. Both questions have been extensively studied for the nonlocal-interaction equations (NLI) which are Wasserstein gradient flows of (1.5) with $P \equiv 0$. For attractive interaction potentials it was shown that the solutions converge to a delta mass [7], while for more general repulsive-attractive potentials very rich families of steady states were discovered [3,35]. We remark that the dynamics of the (NL²IE) can be significantly different. Namely, as the example of Remark 3.18 shows, the solutions for attractive potentials do not necessarily converge to a point.

A further question closely related to asymptotics is the contractivity of solutions of (NL²IE). For Riemannian gradient flows the contractivity of the flow follows from the geodesic convexity of the energy. In particular if $K(x, y) = k(x - y)$, where $k$ is symmetric and convex, the NLI flow is contractive in Wasserstein metric [2,11]. Determining the geodesic convexity of energies in the setting of the nonlocal Wasserstein metrics is an intriguing question. Thus far, the only result in the general (not purely discrete) setting is the geodesic convexity of the entropy [23]. However, for Finslerian gradient flows we caution that establishing geodesic convexity does
not imply contractivity, as [42] shows. Instead a new property of skew-convexity [42, Definition 3.1] needs to be investigated.

Finally we note that the asymptotics of gradient flows with respect to (nonlocal) Wasserstein metrics in discrete setting has recently been investigated in [15, 26], where the equations also include diffusion (i.e., energy includes an entropic contribution). These papers use the convexity of the total energy in the discrete setting to establish the exponential convergence of the flow towards the unique minimizer. Establishing under which conditions (on the graph construction, etc.) do these estimates persist in the discrete to continuum limit as the number of vertices increases is an interesting open problem. We also remark that, while these results do not carry over to our setting, analyzing the asymptotics of (NL²IE) in purely discrete setting is an intriguing and potentially approachable question.

### 1.3. Relation to the Numerical Finite-Volume Upwind Scheme

Equation (1.7) can be interpreted in several ways. For example, it can be understood as the master equation of a continuous-time and continuous-space Markov jump process on the graphon $(\mathbb{R}^d, \eta)$, that is, a continuous graph with vertices $\mathbb{R}^d$, and symmetric weight $\eta(x, y)$ for $(x, y) \in \{(x, y) \in \mathbb{R}^d \times \mathbb{R}^d : x \neq y\}$. The stochastic interpretation is that a particle at position $x \in \mathbb{R}^d$ jumps according to the measure $v(x, y)\eta(x, y) d\mu(y)$ to $y \in \mathbb{R}^d$. In this way it gives rise to a Markov jump process related to the numerical upwind scheme.

The numerical upwind scheme is one of the basic finite-volume methods used to solve conservation laws; see [29]. To draw the connection, let $\{x_1, \dots, x_n\}$ be a suitable representative of a tessellation $\{K_1, \dots, K_n\}$, for instance a Voronoi tessellation, of some bounded domain $\Omega \subset \mathbb{R}^d$. Let $\mu$ be the Lebesgue measure on $\Omega$ and take $\eta$ to be the transmission coefficient common in finite-volume schemes: $\eta(x_i, x_j) = \mathcal{H}^{d-1}(\overline{K_i} \cap \overline{K_j}) / \text{Leb}(K_i)$, for $i, j \in \{1, \dots, n\}$, where $\mathcal{H}^{d-1}(\overline{K_i} \cap \overline{K_j})$ is the $d-1$ dimensional Hausdorff measure of the common face between $K_i$ and $K_j$. With this choice the equation (1.7) becomes the (continuous-time) discretization of the classical continuity equation

$$ \partial_t \rho_t + \nabla \cdot (v_t \rho_t) = 0 $$

for some vector field $v_t: \Omega \to \mathbb{R}^d$. Hereby, the discretized vector field $v_t$ is obtained from $v_t$ by taking the average over common interfaces:

$$ v_t(x_i, x_j) = \frac{1}{\mathcal{H}^{d-1}(\overline{K_i} \cap \overline{K_j})} \int_{\overline{K_i} \cap \overline{K_j}} v_t \cdot v_{K_i, K_j} \, d\mathcal{H}^{d-1}, $$

where $v_{K_i, K_j}$ is the unit normal to $K_i$ pointing from $K_i$ to $K_j$. We refer to the recent work [9] for a variational interpretation of the upwind scheme, which is close to that we propose for the more general equation (1.7). Earlier results in this direction are contained in [21, 38].

The connection to finite-volume schemes explains also that the nonlocality in (1.7) introduces a regularization, which in the numerical literature is referred to as numerical diffusion. That the numerical diffusion is actually an honest Markov jump process, as described at the beginning of this section, was observed and used to find optimal convergence rates in the works [19, 20, 45, 46].
## 1.4. Comparison with Other Discrete Metrics and Gradient Structures

The interpretation of diffusion on graphs as gradient flows of the entropy was independently carried out in [14, 36, 37]. Here we recall the descriptions of the flows relying on reversible Markov chains, which was the framework used in [25, 27, 36]. Starting with Markov chains, which then determine the edge weights, offers an additional layer of modeling flexibility. In particular, consider the Markov chain with state space $X = \{x_1, \dots, x_n\}$ and jump rates $\{Q_{x,y}\}_{x,y \in X}$. Let $\pi_x$ be the reversible probability measure for the Markov chain, meaning that it satisfies the detailed balance condition $\pi_x Q_{x,y} = \pi_y Q_{y,x}$. The edge weights $\{w_{x,y}\}_{x,y \in X}$ are given by $w_{x,y} = \pi_x Q_{x,y}$. The energy considered is the relative entropy: for $\rho: X \to [0, 1]$ with $\sum_{x \in X} \rho_x = 1$ we define

$$ \mathcal{H}(\rho | \pi) = \sum_x \rho_x \log \frac{\rho_x}{\pi_x} = \sum_x \rho_x \log \rho_x - \sum_x \rho_x \log \pi_x = \mathcal{S}(\rho) + \mathcal{E}_P(\rho) \tag{1.9} $$
with $P_x = -\log \pi_x$.

The paths in the configuration space are given as the solution of the continuity equation which for the flux $\{j_{x,y} : [0, T] \to \mathbb{R}\}_{x,y \in X}$ takes the form (1.2).

To compute the flux from a given velocity $\{v_{x,y}\}_{x,y \in X}$ (an edge-based quantity) and density $\{\rho_x\}_{x \in X}$ (a vertex-based quantity), one interpolates the densities at vertices to define the density (and hence the flux) along the edges. The literature so far has considered a proportional constitutive relation of the form

$$ j_{x,y} = v_{x,y} \theta \left( \frac{\rho_x}{\pi_x}, \frac{\rho_y}{\pi_y} \right), \tag{1.10} $$

where the function $\theta: \mathbb{R}_+ \times \mathbb{R}_+ \to \mathbb{R}_+$ needs to be one-homogeneous for dimensional reasons. In addition, it is assumed that the function $\theta$ is an interpolation, that is, $\min\{a, b\} \le \theta(a, b) \le \max\{a, b\}$. The choice providing a gradient flow characterization for linear Markov chains is the logarithmic mean, defined by $\theta(a, b) = \frac{a-b}{\log a - \log b}$ for $a \ne b$ and $\theta(a, a) = a$.

The associated transportation distance is obtained by minimizing the action functional

$$ \mathcal{A}(\rho, j) = \frac{1}{2} \sum_{x,y} j_{x,y} v_{x,y} w_{x,y} = \frac{1}{2} \sum_{x,y} \frac{|j_{x,y}|^2}{\theta\left(\frac{\rho_x}{\pi_x}, \frac{\rho_y}{\pi_y}\right)} \pi_x Q_{x,y}. \tag{1.11} $$

The corresponding transportation distance is induced as the minimum of the action along paths:

$$ \inf \left\{ \int_0^1 \mathcal{A}(\rho(t), j(t)) dt : (\rho(t), j(t))_{t \in [0,1]} \text{ solves (1.2) and } \rho(0) = \rho_0, \rho(1) = \rho_1 \right\}. $$

As we do in Corollary 2.8, it was shown that it suffices to consider antisymmetric fluxes. To arrive at a gradient flow formulation, one considers the metric induced by the action function (1.11):
$$g_{\rho}(j^1, j^2) = \frac{1}{2} \sum_{x,y} \frac{j_{x,y}^1 j_{x,y}^2}{\theta\left(\frac{\rho_x}{\pi_x}, \frac{\rho_y}{\pi_y}\right)} \pi_x Q_{x,y}. \quad (1.12)$$

Then the gradient grad $\mathcal{H}$ of the relative entropy (1.9) with respect to this metric is given as the antisymmetric flux $j^*$ of minimal norm satisfying

$$g_{\rho}(\text{grad } \mathcal{H}, j) = \text{Diff } \mathcal{H}[j] = \left. \frac{d}{dt} \right|_{t=0} \mathcal{H}(\tilde{\rho}(t)), \quad (1.13)$$

for any curve $(\tilde{\rho}(t))_{t \ge 0}$ such that $\partial_t \rho(0) = -(\bar{\nabla} \cdot j)$. Expanding (1.13) and using that $j^*$ is antisymmetric gives

$$\frac{1}{2} \sum_{x,y} \frac{j_{x,y} j_{x,y}^*}{\theta\left(\frac{\rho_x(t)}{\pi_x}, \frac{\rho_y(t)}{\pi_y}\right)} \pi_x Q_{x,y} = -\frac{1}{2} \sum_{x,y} \left( \log \frac{\rho_x(t)}{\pi_x} - \log \frac{\rho_y(t)}{\pi_y} \right) j_{x,y} w_{x,y}.$$

Since this identity holds for all $j_{x,y}$, the flux $j^*$ is identified by

$$j_{x,y}^* = - \left( \log \frac{\rho_x}{\pi_x} - \log \frac{\rho_y}{\pi_y} \right) \theta \left( \frac{\rho_x(t)}{\pi_x}, \frac{\rho_y(t)}{\pi_y} \right) = - \left( \frac{\rho_x(t)}{\pi_x} - \frac{\rho_y(t)}{\pi_y} \right),$$

where the last equality holds for the particular choice of the logarithmic mean interpolation $\theta(r, s) = \frac{r-s}{\ln r - \ln s}$. By plugging $j_{x,y}^*$ into the continuity equation (1.2), one recovers the (linear) heat equation on graphs.

The next relevant step is the introduction of the interaction and the potential energies as in (1.1). In particular, [25] provides a gradient flow structure for free energy functionals of the form

$$\mathcal{F}_{\beta}(\rho) = \beta^{-1} S(\rho) + \mathcal{E}_{X}(\rho), \quad (1.14)$$

where $\beta > 0$ is the inverse temperature. This is the discrete analogue of the McKean-Vlassov equation. Finding a desirable gradient flow structure is nontrivial since considering the logarithmic interpolation, which makes the diffusion term linear, would make the potential term nonlinear, and thus the Fokker-Planck equation on graphs would be nonlinear. To cope with this, the framework of [25] extends the linear theory outlined above to a family of nonlinear Markov chains satisfying a local detailed balance condition. The consequence for the resulting gradient structure is that the quantities $\{\pi_x\}_{x \in X}$, $\{Q_{x,y}\}_{x,y \in X}$ and $\{w_{x,y}\}_{x,y \in X}$ depend on the current state $\rho$ in such a way that the detailed balance condition $w_{x,y}[\rho] = \pi_x[\rho]Q_{x,y}[\rho] = \pi_y[\rho]Q_{y,x}[\rho]$ is still valid for all $\rho \in \mathcal{P}(X)$. In particular, for $\mathcal{F}_{\beta}$ defined in (1.14), it holds that

$$\pi_x[\rho] = \frac{1}{Z_{\beta}} \exp \left( -\beta \left( P_x + \sum_y K_{x,y} \rho_y \right) \right) \quad \text{with} \\ Z_{\beta} = \sum_x \exp \left( -\beta \left( P_x + \sum_y K_{x,y} \rho_y \right) \right).$$
It would be natural to try to build the framework for the case $\beta = \infty$, which we consider in this paper, by taking the limit $\beta \to \infty$ in the framework of [25]. It turns out that this limit is singular for the constructed gradient structure. First of all, the measure $\pi_x[\rho]$ degenerates at all points except at the argmin of the effective potential $x \mapsto P_x + \sum_y K_{x,y}\rho_y$. This causes the constitutive relation (1.10) to become meaningless. A more detailed analysis also shows that the metric in (1.12) degenerates.

We also note that in this setting the potential functions $P$ and $K$ and inverse temperature $\beta$ enter the metric in (1.11) through the weights $w_{x,y}$ and rate matrix $Q_{x,y}$. This is in stark contrast to the continuous classical gradient flow formulation for free energies of the form $\mathcal{F}_\beta$ form (1.14), where the metric is always the $L^2$-Wasserstein distance, independently of the potentials $P$ and $K$ and also of the inverse temperature $\beta > 0$, including $\beta = \infty$ [2, 10, 11, 33].

Another approach to McKean-Vlasov equations is to consider the arithmetic interpolation, as was done in [15]. The theory the authors developed requires the densities to be strictly positive and diffusion to be present. We note that the diffusion itself is nonlinear.

The above problems lead us to consider the upwind interpolation in the flux-velocity relation (1.10). In view of (1.2), this relation is replaced in the present setting by

$$j_{x,y} = \rho_x(v_{x,y})_+ - \rho_y(v_{x,y})_- = \Theta(\rho_x, \rho_y; v_{x,y})v_{x,y} \quad \text{where}$$

$$\Theta(a, b; v) = \begin{cases} a & \text{if } v > 0, \\ b & \text{if } v < 0, \\ 0 & \text{if } v = 0. \end{cases} \qquad (1.15)$$

Note that the relation (1.15) is a functional relation between velocity and flux with the interpolation $\Theta$ depending on the velocity.

We remark that solutions of system (1.2)–(1.4) are not the limit of the gradient flows in [25] as $\beta \to \infty$. We emphasize here that the limit of these dynamics as $\beta \to \infty$ would in fact not be the desirable gradient flow of the nonlocal-interaction energy, since the initial support of the solutions would never expand; see the related Remark 1.2.

We conclude this section by observing that it seems possible to generalize the upwind interpolation in a continuous way to define a flux-velocity relation to deal with free energies $\mathcal{F}_\beta$ for $\beta > 0$. A candidate, inspired by the Scharfetter–Gummel scheme [44], is the following constitutive flux-velocity relation depending on $\beta$:

$$j_{x,y}^{\beta} = v_{x,y} \frac{\rho_x \exp(\beta v_{x,y}/2) - \rho_y \exp(-\beta v_{x,y}/2)}{\exp(\beta v_{x,y}/2) - \exp(-\beta v_{x,y}/2)}.$$

In particular, it holds that $j_{x,y}^{\beta} \to j_{x,y}$ as $\beta \to \infty$, where $j_{x,y}$ is as in (1.15). The form of $j_{x,y}^{\beta}$ can be physically deduced from the one-dimensional cell problem for the unknown value $j_{x,y}^{\beta} \in \mathbb{R}$ and function $\rho: [0, 1] \to \mathbb{R}$:

$$j_{x,y}^{\beta} = -\beta^{-1} \bar{\nabla} \rho(\cdot) + v_{x,y} \rho(\cdot) \quad \text{on } [0, 1], \quad \text{with } \rho(0) = \rho_x \text{ and } \rho(1) = \rho_y.$$
Note that $j_{x,y}^{\beta} = \frac{\rho_x - \rho_y}{\beta}$ for $v_{x,y} = 0$, which is the flux due to Fick's law. Likewise, $j_{x,y}^{\beta} = 0$ for $v_{x,y} = \beta^{-1} \log \frac{\rho_y}{\rho_x}$, which is the velocity needed to counteract the diffusion. In [47], it is shown that the Scharfetter-Gummel finite volume scheme provides a stable positivity preserving numerical approximation of the diffusion-aggregation equation, which also respects the thermodynamic free energy structure. We pursue the investigation of the existence of a possible related gradient structure in future research.

## 1.5. Connections to Machine Learning

Part of the motivation for the present work comes from applications to machine learning. Here we introduce a family of nonlinear gradient flows that is relevant to discovering local concentrations in networks akin to modes of a distribution.

Our main interest is in equations posed on graphs whose vertices are random samples of some underlying distribution and whose edge weights are a function of distances between vertices. In machine learning one often deals with data in the form of a point cloud in high-dimensional space. While the ambient dimension may be very large, the data often possess an underlying low-dimensional structure that can be used in making reliable inferences about the underlying data distribution. To use the geometric information, we follow one of the standard approaches and consider graphs associated to point clouds. Formulating the machine learning tasks directly on the point cloud enables one to access the geometric structure of the distribution in a simple and computationally efficient way. The works in the literature have mostly focused on models based on minimizing objective functionals modeling tasks such as clustering or dimensionality reduction [5,31,32,34,40], or based on characterizing clusters through estimating some property of the data distribution (most often the density); see [12] and references therein. Only few dynamical models have been considered—notable among them are diffusion maps [16], where the heat equation is used to redistance the points.

Here we focus on models that are motivated by nonlocal PDEs. Consider a probability measure $\mu$ on $\mathbb{R}^d$ with finite second moments. Let $X = \{x_1, \dots, x_n\}$ be random i.i.d. samples of the measure $\mu$. Let $\mu^n = \frac{1}{n} \sum_{i=1}^n \delta_{x_i}$ be the empirical measure of the sample and let $K: \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}$ be symmetric and $P: \mathbb{R}^d \to \mathbb{R}$. The total energy $E_X: P(X) \to \mathbb{R}$, given in (1.1), for the empirical measure $\mu^n$ can be rewritten as

$$ \begin{aligned} E_X(\mu^n) &= E_I(\mu^n) + E_P(\mu^n) \\ &= \frac{1}{2} \iint_{\mathbb{R}^d \times \mathbb{R}^d} K(x, y) \, d\mu^n(x) \, d\mu^n(y) + \int_{\mathbb{R}^d} P(x) \, d\mu^n(x). \end{aligned} \quad (1.16) $$

The gradient flow of $E_X$ with respect to the graph Wasserstein metric $T_{\mu^n}$ defined in (1.8) is described by the ODE system (1.2)–(1.4), where $K_{x_i, x_j} = K(x_i, x_j)$ and $P_{x_i} = P(x_i)$ for all $i, j \in \{1, \dots, n\}$. Another evolution by such system is illustrated on Fig. 2.

Here we remark on the contrast between (1.2)–(1.4) and the gradient flow of (1.16) in the ambient space $\mathbb{R}^d$, with respect to the standard Wasserstein metric,
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** A solution of the nonlocal-interaction equation on graphs driven by the energy (1.1). We consider a random geometric graph based on 240 sample points X from a 2D bean data set. The connectivity distance is $\varepsilon = 0.23$. The edge weights are $w_{x,y} = \exp(-24|x-y|^2)$, provided that the vertices x and y of X are connected. The interaction potential is $K_{x,y} = 1 - \exp(-8|x-y|^2)$ and the external potential is $P \equiv 0$. The solution, starting from a uniform distribution, is shown at time $t = 200$. Brighter color indicates more mass (right)</figcaption></figure>

which takes the form

$$ \dot{x}_i = -\nabla P(x_i) - \sum_{j=1}^{n} \rho_j \nabla_{x_i} K(x_i, x_j). \quad (1.17) $$

The first notable difference is that, on the graph, masses change and the positions remain fixed, while in $\mathbb{R}^d$ positions change and the masses remain fixed. This difference is somewhat superficial, since both equations describe the rearrangement of mass in order to decrease the same energy in the most efficient way measured by two different metrics. The main difference is that the graph encodes the geometry of the space that mass is allowed to occupy. In particular, it ensures that the geometric mode discovered will be a data point itself.

We note that the popular mean-shift algorithm [17] can be interpreted as a time-stepping algorithm to approximate solutions of (1.17) with $K \equiv 0$ and $P = \ln(\theta * \mu^n(0))$, where $\mu^n(0)$ is the empirical measure of the initial distribution of particles and $\theta * \mu^n(0)$ is the kernel density estimate of the density $\rho$ of the underlying distribution. Namely the step of the mean-shift algorithm is to replace the position of the particle at $x_j$ by the center of mass of $\theta(\cdot - x_j) * \mu_n(0)$ and iterate the procedure. Formal expansion shows that this is a time step of the forward scheme for the flow driven by $P = \ln(\theta * \mu^n(0))$. We note that considering the gradient flow of the corresponding energy on the graph (1.2)–(1.4) ensures that the modes of the distribution discovered by the (graph) mean-shift algorithm will remain within the data set. Furthermore, we note that adding nonlocal attraction on the graph progressively clumps nearby masses together and thus provides an approach to agglomerative clustering.

One of our main results, stated in Theorem 3.14, is that as $n \to \infty$ the solutions of the graph-based equation (1.2)–(1.4) narrowly converge along a subsequence to a solution of the nonlocal nonlocal-interaction equation (NL$^2$IE).
## 2. Nonlocal Continuity Equation and Upwind Transportation Metric

### 2.1. Weight Function

Throughout the paper we consider a weight function $\eta: \{(x, y) \in \mathbb{R}^d \times \mathbb{R}^d : x \neq y\} \to [0, \infty)$, which shall always satisfy

$$ 
\left\{ 
\begin{array}{l} 
\eta \text{ is continuous on } G = \{(x, y) \in \mathbb{R}^d \times \mathbb{R}^d : x \neq y, \eta(x, y) > 0\}; \\ 
\eta \text{ is symmetric and non-negative, that is } \eta(x, y) = \eta(y, x) \text{ for all } (x, y) \in G. 
\end{array} 
\right. 
\qquad (W) 
$$

Since $\eta$ is symmetric, we regard the edges set $G$ as undirected graph. Many of the edge-based quantities we consider, like vector fields and fluxes, will lie in an $\eta$-weighted $L^2$ space, $L^2(\eta \lambda)$ for some $\lambda \in \mathcal{M}(G)$. The space $L^2(\eta \lambda)$ is equipped with the inner product

$$ \langle f, g \rangle_{L^2(\eta\lambda)} = \frac{1}{2} \iint_G f(x, y)g(x, y)\eta(x, y)d\lambda(x, y) \quad \text{for all } f, g \in L^2(\eta\lambda), \quad (2.1) $$

where the factor $\frac{1}{2}$ ensures that each undirected edge is counted only once.

Below we state two assumptions on the base measure $\mu \in \mathcal{M}^+(\mathbb{R}^d)$ and the weight function $\eta$, where we use the notation $\vee$ to denote the maximum.

(A1) (moment bound) The family of functions $\{|x - \cdot|^2 \vee |x - \cdot|^4\} \eta(x, \cdot)_{x \in \mathbb{R}^d}$ is uniformly integrable with respect to $\mu$, that is, for some $C_\eta \in (0, \infty)$, it holds that

$$ \sup_{x \in \mathbb{R}^d} \int \left( |x - y|^2 \vee |x - y|^4 \right) \eta(x, y) \, d\mu(y) \leq C_\eta. $$

(A2) (local blow-up control) The family of measures $\{|x - \cdot|^2 \eta(x, \cdot) \mu(\cdot)\}_{x \in \mathbb{R}^d}$ is locally uniformly integrable, that is,

$$ \lim_{\epsilon \to 0} \sup_{x \in \mathbb{R}^d} \int_{B_{\epsilon}(x) \setminus \{x\}} |x - y|^2 \eta(x, y) \, d\mu(y) = 0, \quad \text{where} \\ B_{\epsilon}(x) = \{y \in \mathbb{R}^d : |x - y| < \epsilon\}. $$

**Remark 2.1.** Continuity on $G$ in (W) is needed to obtain lower semicontinuity of the action functional; see Lemma 2.9. Assumption (A1) ensures well-posedness of the nonlocal continuity equation we shall introduce in Section 2.3, whereas Assumption (A2) is necessary for compactness of solutions to the nonlocal continuity equation; see Proposition 2.17.

**Example 2.2.** Typically the function $\eta$ is a function of the distance

$$ \eta(x, y) = \vartheta(|x - y|) \quad \text{for all } (x, y) \in G, $$

where $\vartheta: (0, \infty) \to [0, \infty)$ is continuous on $\{\vartheta > 0\}$ and satisfies analogues of (A1) and (A2). An important example are geometric graphs with connectivity distance given by $\epsilon > 0$ and weight

$$ \eta_{\epsilon}(x, y) = \frac{2(2 + d) \chi_{B_{\epsilon}(x)}(y)}{\epsilon^2 |B_{\epsilon}|} \quad \text{for all } (x, y) \in G. \quad (2.2) $$
In this example, fixing $\mu = \text{Leb}(\mathbb{R}^d)$, we conjecture that the weak formulation of (NL$^2$IE)—see Section 3—converges to the nonlocal aggregation equation $\partial_t \rho_t = \nabla \cdot (\rho_t \nabla K * \rho_t + \rho_t \nabla P)$ as $\varepsilon \to 0$ for sufficiently smooth potentials K and P. See Section 3.5 for a discussion on the local limit.

## 2.2. Action

The form of the action inside (1.8) seems practical, but it does not have any obvious convexity and lower semicontinuity properties. Therefore, we define the action in flux variables. We start by introducing some notation. For a signed measure $j \in \mathcal{M}(G)$, we denote by $j = j^+ - j^-$ its Jordan decomposition. Moreover, for any measurable $A \subseteq G$, let $A^\top = \{(y, x) \in \mathbb{R}^d \times \mathbb{R}^d : (x, y) \in A\}$ be its transpose. Likewise, for $j \in \mathcal{M}(G)$, we denote by $j^\top$ the transposed measure defined by $j^\top(A) = j(A^\top)$.

For any measures $\mu \in \mathcal{M}^+(\mathbb{R}^d)$ and $\rho \in \mathcal{P}(\mathbb{R}^d)$, we define the (restricted) product measures $\gamma_i \in \mathcal{M}^+(G)$ for $i = 1, 2$ as

$$ \begin{aligned} & d\gamma_1(x, y) = d\rho(x) d\mu(y) \quad \text{and} \\ & d\gamma_2(x, y) = d\mu(x) d\rho(y) \quad \text{for } (x, y) \in G. \end{aligned} \quad (2.3) $$

Note that $\gamma_1^\top = \gamma_2$. We define the action for general $\eta$ which we only require to satisfy Assumption (W), i.e., continuity on $G$, symmetry and positivity.

**Definition 2.3. (Action)** For $\mu \in \mathcal{M}^+(\mathbb{R}^d)$, $\rho \in \mathcal{P}(\mathbb{R}^d)$ and $j \in \mathcal{M}(G)$, consider $\lambda \in \mathcal{M}(G)$ such that $\rho \otimes \mu, \mu \otimes \rho, |j| \ll |\lambda|$. We define

$$ \begin{aligned} \mathcal{A}(\mu; \rho, j) = & \frac{1}{2} \iint_G \alpha \left( \frac{d\mathbf{j}}{d|\lambda|}, \frac{d(\rho \otimes \mu)}{d|\lambda|} \right) \eta \, d|\lambda| \\ & + \frac{1}{2} \iint_G \alpha \left( -\frac{d\mathbf{j}}{d|\lambda|}, \frac{d(\mu \otimes \rho)}{d|\lambda|} \right) \eta \, d|\lambda|. \end{aligned} \quad (2.4) $$

Hereby, the lower semicontinuous, convex, and positively one-homogeneous function $\alpha: \mathbb{R} \times \mathbb{R}_+ \to \mathbb{R}_+ \cup \{\infty\}$ is defined, for all $j \in \mathbb{R}$ and $r \ge 0$, by

$$ \alpha(j, r) := \begin{cases} \frac{(j_+)^2}{r} & \text{if } r > 0, \\ 0 & \text{if } j \le 0 \text{ and } r = 0, \\ \infty & \text{if } j > 0 \text{ and } r = 0, \end{cases} \quad (2.5) $$

with $j_+ = \max\{0, j\}$. If the measure $\mu$ is clear from the context, we write $\mathcal{A}(\rho, j)$ for $\mathcal{A}(\mu; \rho, j)$.

Note that Definition 2.3 is well-posed since the one-homogeneity of $\alpha$ makes it independent of the particular choice of $\lambda$ as long as the absolute continuity condition in Definition 2.3 is satisfied. An example of such measure is a $\lambda$ such that $|\lambda| = |\rho \otimes \mu| + |\mu \otimes \rho| + |j|$. Moreover, $\lambda$ can be chosen symmetric, otherwise it can be replaced by $\frac{1}{2}(\lambda + \lambda^\top)$.
**Remark 2.4.** We note that the action is inversely proportional to the measure $\mu$: doubling the measure $\mu$ leads to halving the action. This has important consequence for the way $\mu$ influences the geometry of the space of measures. In particular, $\mu$ not only sets the region where mass can be transported, but also makes the transport less costly in the regions of high density of $\mu$.

**Remark 2.5.** If $\rho \ll \mu$, then we denote its density by $\rho$ by abuse of notation, and if furthermore $j \ll \mu \otimes \mu$ with density $j$, then it holds that

$$A(\mu; \rho, j) = \frac{1}{2} \iint_G \left( \frac{(j(x, y))_+^2}{\rho(x)} + \frac{(j(x, y))_-^2}{\rho(y)} \right) \eta(x, y) \, d\mu(x) \, d\mu(y).$$

In the following lemma we can see that the action takes the form from the tentative definition of the metric in (1.8), as soon as it is bounded.

**Lemma 2.6.** Let $\mu \in \mathcal{M}^+(\mathbb{R}^d)$, $\rho \in \mathcal{P}(\mathbb{R}^d)$ and $j \in \mathcal{M}(G)$ be such that $A(\mu; \rho, j) < \infty$. Then there exists a measurable $v: G \to \mathbb{R}$ such that

$$d\boldsymbol{j}(x, y) = v(x, y)_+ d\rho(x) d\mu(y) - v(x, y)_- d\mu(x) d\rho(y), \quad (2.6)$$

and it holds that

$$A(\mu; \rho, j) = \frac{1}{2} \iint_G (|v(x, y)_+|^2 + |v(x, y)_-|^2) \eta(x, y) \, d\rho(x) \, d\mu(y). \quad (2.7)$$

In particular, if $v \in \mathcal{V}^{\text{as}}(G)$, then

$$A(\mu; \rho, j) = \iint_G |v(x, y)_+|^2 \eta(x, y) \, d\rho(x) \, d\mu(y). \quad (2.8)$$

**Proof.** Let $\lambda \in \mathcal{M}^+(G)$ be such that $d\gamma_1(x, y) = d\rho(x) d\mu(y) = \tilde{\gamma}_1(x, y) d\lambda(x, y)$, likewise $d\gamma_2(x, y) = d\mu(x) d\rho(y) = \tilde{\gamma}_2(x, y) d\lambda(x, y)$, and $d\boldsymbol{j} = \tilde{j} d\lambda$ for some measurable $\tilde{\gamma}_1, \tilde{\gamma}_2, \tilde{j}: G \to \mathbb{R}$. Without loss of generality we can assume $\lambda$ to be symmetric; for instance by considering $\frac{1}{2}(\lambda + \lambda^\top)$ instead. Thus, (2.4) implies

$$A(\mu; \rho, j) = \frac{1}{2} \iint_G \left( \alpha(\tilde{j}(x, y), \tilde{\gamma}_1(x, y)) + \alpha(-\tilde{j}(x, y), \tilde{\gamma}_2(x, y)) \right) \eta(x, y) d\lambda(x, y) < \infty.$$

By the definition of the function $\alpha$ in (2.5), it immediately follows that the vector field $\tilde{v}^+(x, y) = \frac{\tilde{j}(x, y)_+}{\tilde{\gamma}_1(x, y)}$ is well-defined $\gamma_1$-a.e. on $G$. By the same argument, we find that $\tilde{v}^-(x, y) = \frac{\tilde{j}(x, y)_-}{\tilde{\gamma}_2(x, y)}$ is well-defined $\gamma_2$-a.e. on $G$. Since $\gamma_1 = \gamma_2^\top$ we have that $(\tilde{v}^-)^\top$ exists $\gamma_1$-a.e. on $G$. Hence, we obtain the measurable vector field

$$v(x, y) = \tilde{v}^+(x, y) - \tilde{v}^-(x, y).$$

The statement (2.7) follows by using the positively one-homogeneity of $\alpha$, the identity $\alpha(j, r) = \alpha(j_+, r)$ and the symmetry of $\lambda$:

$$A(\mu; \rho, j) = \frac{1}{2} \iint_G \alpha(v(x, y)_+ \tilde{\gamma}_1(x, y), \tilde{\gamma}_1(x, y)) \eta(x, y) d\lambda(x, y)$$
$$ 
\begin{aligned}
&+\frac{1}{2}\iint_G \alpha(v(x, y)-\tilde{\gamma}_2(x, y), \tilde{\gamma}_2(x, y)) \eta(x, y) \, d\lambda(x, y) \\
&=\frac{1}{2}\iint_G |v(x, y)|^2 \tilde{\gamma}_1(x, y) \eta(x, y) \, d\lambda(x, y) \\
&+\frac{1}{2}\iint_G |v(y, x)|^2 \tilde{\gamma}_1(x, y) \eta(x, y) \, d\lambda(x, y).
\end{aligned}
\quad \square
$$

**Definition 2.7.** (Nonlocal gradient and divergence) For any function $\phi: \mathbb{R}^d \to \mathbb{R}$ we define its nonlocal gradient $\bar{\nabla}\phi: G \to \mathbb{R}$ by

$$ \bar{\nabla}\phi(x, y) = \phi(y) - \phi(x) \quad \text{for all } (x, y) \in G. $$

For any $j \in \mathcal{M}(G)$, its nonlocal divergence $\bar{\nabla} \cdot j \in \mathcal{M}(\mathbb{R}^d)$ is defined as $\eta$-weighted adjoint of $\bar{\nabla}$, i.e.,

$$ 
\begin{aligned}
\int \phi \, d\bar{\nabla} \cdot j &= -\frac{1}{2} \iint_G \bar{\nabla} \phi(x, y) \eta(x, y) \, dj(x, y) \\
&= \frac{1}{2} \int \phi(x) \int \eta(x, y) (dj(x, y) - dj(y, x)).
\end{aligned}
$$

In particular, for $j \in \mathcal{M}^{\text{as}}(G) := \{j \in \mathcal{M}(G) : j^\top = -j\}$,

$$ \int \phi \, d\bar{\nabla} \cdot j = \iint_G \phi(x) \eta(x, y) \, dj(x, y). $$

If $j$ is given by (2.6) for some $v \in \mathcal{V}^{\text{as}}(G)$, then the flux satisfies an antisymmetric relation on the support of $\gamma_1$-a.e. on $G$, i.e., $j^+ = (j^\top)^{-} \gamma_1$-a.e. on $G$. The following corollary shows that those antisymmetric fluxes are the relevant ones for the minimization of the action functional. For this reason, the natural class of fluxes are those measure on $G$ which are antisymmetric with positive part absolutely continuous with respect to $\gamma_1$, that is,

$$ \mathcal{M}_{\gamma_1}^{\text{as}}(G) = \{j \in \mathcal{M}(G) : j^+ \ll \gamma_1, j^- \ll \gamma_1^\top, j^+ = (j^\top)^{-} \gamma_1 \text{-a.e.}\} \quad (2.9) $$

**Corollary 2.8.** (Antisymmetric vector fields have lower action) Let $\mu \in \mathcal{M}^+(\mathbb{R}^d)$, $\rho \in \mathcal{P}(\mathbb{R}^d)$ and $j \in \mathcal{M}(G)$ be such that $\mathcal{A}(\mu; \rho, j) < \infty$. Then there exists an antisymmetric flux $j^{\text{as}} \in \mathcal{M}_{\gamma_1}^{\text{as}}$ such that

$$ \bar{\nabla} \cdot j = \bar{\nabla} \cdot j^{\text{as}}, $$

with lower action:

$$ \mathcal{A}(\mu; \rho, j^{\text{as}}) \le \mathcal{A}(\mu; \rho, j). $$
**Proof.** Let us set $j^{\text{as}} = (j - j^{\top})/2$. Since $\eta$ is symmetric and $(\bar{\nabla}\phi)^{\top} = -\bar{\nabla}\phi$, we get

$$ \begin{aligned} \iint_G \bar{\nabla}\phi \eta \,d\mathbf{j}^{\text{as}} &= \frac{1}{2}\iint_G \bar{\nabla}\phi \eta (d\mathbf{j} - d\mathbf{j}^{\top}) \\ &= \frac{1}{2}\iint_G \bar{\nabla}\phi \eta \,d\mathbf{j} - \frac{1}{2}\iint_G (\bar{\nabla}\phi)^{\top} \eta \,d\mathbf{j} \\ &= \iint_G \bar{\nabla}\phi \eta \,d\mathbf{j}. \end{aligned} $$

By an application of Lemma 2.6 and comparison of (2.7) and (2.8) it is enough to show that, for all $(x, y) \in G$,

$$ \begin{aligned} &|v^{\text{as}}(x, y)_+|^2 + |v^{\text{as}}(x, y)_-|^2 + |v^{\text{as}}(y, x)_+|^2 + |v^{\text{as}}(y, x)_-|^2 \\ &\le |v(x, y)_+|^2 + |v(x, y)_-|^2 + |v(y, x)_+|^2 + |v(y, x)_-|^2 \end{aligned} $$

for any measurable $v: G \to \mathbb{R}$, where $v^{\text{as}}(x, y) = (v(x, y) - v(y, x))/2$. This estimate is a consequence of Jensen's inequality applied to the convex functions

$$ \varphi^{\pm}: \mathbb{R} \to \mathbb{R} \quad \text{with} \quad \varphi^{\pm}(r) = (r_{\pm})^2. \quad \square $$

**Lemma 2.9.** (Lower semicontinuity of the action) The action is lower semicontinuous with respect to the narrow convergence in $M^+(\mathbb{R}^d) \times \mathcal{P}(\mathbb{R}^d) \times \mathcal{M}(G)$. That is, if $\mu^n \to \mu$ in $\mathcal{M}(\mathbb{R}^d)$, $\rho^n \to \rho$ in $\mathcal{P}(\mathbb{R}^d)$, and $j^n \to j$ in $\mathcal{M}(G)$, then

$$ \liminf_{n \to \infty} \mathcal{A}(\mu^n; \rho^n, j^n) \ge \mathcal{A}(\mu; \rho, j). $$

**Proof.** First, note that the narrow convergence of any sequences $(\rho^n)_n$ and $(\mu^n)_n$ implies the narrow convergence of the product: $\rho^n \otimes \mu^n \to \rho \otimes \mu$ in $\mathcal{P}(\mathbb{R}^d) \times M^+(\mathbb{R}^d)$, therefore also in $M^+(G)$. Then, in Definition 2.3 consider the vector-valued measure

$$ \lambda = (j, \rho \otimes \mu, \mu \otimes \rho). $$

Further, we define the function

$$ f: G \times \mathbb{R}^3 \to \mathbb{R} \quad \text{by} \quad f((x, y), (j, \gamma_1, \gamma_2)) = (\alpha(j, \gamma_1) + \alpha(-j, \gamma_2)) \eta(x, y). $$

Since the function $\eta$ is lower semicontinuous by (W) and $\alpha$ defined in (2.5) is lower semicontinuous, jointly convex and positively one-homogeneous, $f$ satisfies the assumptions of [8, Theorem 3.4.3], whence the claim follows. $\square$

According to Definition 2.3, fluxes and action are strictly related. In case $\mathcal{A}(\mu; \rho, j) < +\infty$, we get a useful upper bound in the following lemma that will be crucial in several technical parts later on.

**Lemma 2.10.** For any $\mu \in M^+(\mathbb{R}^d)$, $\rho \in \mathcal{P}(\mathbb{R}^d)$, $j \in \mathcal{M}(G)$ and any measurable $\Phi: G \to \mathbb{R}_+$ it holds

$$ \left(\frac{1}{2} \iint_G \Phi \eta \,d|\mathbf{j}|\right)^2 \le \mathcal{A}(\mu; \rho, \mathbf{j}) \iint_G \Phi^2 \eta (d\gamma_1 + d\gamma_2). \quad (2.10) $$
**Proof.** Let $\mu \in M^+(\mathbb{R}^d)$, $\rho \in P(\mathbb{R}^d)$ and $j \in M(G)$ be such that $A(\mu; \rho, j) < +\infty$. Let $|\lambda| \in M^+(G)$ be such that $\gamma_1, \gamma_2, |j| \ll |\lambda|$ as in Definition 2.3 and write $\gamma_i = \tilde{\gamma}_i|\lambda|$ and $|j| = |j||\lambda|$ for the densities.

We have that $A := \{(x, y) \in G : \alpha(j, \tilde{\gamma}_1) = \infty \text{ or } \alpha(-j, \tilde{\gamma}_2) = \infty\}$ is a $\lambda$-nullset. We observe the elementary inequality

$$ (j_+)^2 + (j_-)^2 \le \max\{\tilde{\gamma}_1, \tilde{\gamma}_2\}(\alpha(j, \tilde{\gamma}_1) + \alpha(-j, \tilde{\gamma}_2)), \quad \lambda\text{-a.e. in } A^c. $$

In particular, it holds that

$$ |j| = j_+ + j_- \le \sqrt{2 \max\{\tilde{\gamma}_1, \tilde{\gamma}_2\} \sqrt{\alpha(j, \tilde{\gamma}_1) + \alpha(-j, \tilde{\gamma}_2)}}, \quad \lambda\text{-a.e. in } A^c. $$

Hence we can estimate

$$ \begin{aligned} \frac{1}{2} \iint_G \Phi \eta \, d|j| &= \frac{1}{2} \iint_G \Phi \eta |j| \, d|\lambda| = \frac{1}{2} \iint_{A^c} \Phi \eta (j_+ + j_-) \, d|\lambda| \\ &\le \frac{1}{2} \iint_{A^c} \Phi \eta \sqrt{2 \max\{\tilde{\gamma}_1, \tilde{\gamma}_2\}} \sqrt{\alpha(j, \tilde{\gamma}_1) + \alpha(-j, \tilde{\gamma}_2)} \, d|\lambda| \\ &\le \left( \iint_G \Phi^2 \eta \max\{\tilde{\gamma}_1, \tilde{\gamma}_2\} \, d|\lambda| \right)^{\frac{1}{2}} \\ &\quad \times \left( \frac{1}{2} \iint_G (\alpha(j, \tilde{\gamma}_1) + \alpha(-j, \tilde{\gamma}_2)) \eta \, d|\lambda| \right)^{\frac{1}{2}}. \end{aligned} $$

Now, the result follows by estimating $\max\{\tilde{\gamma}_1, \tilde{\gamma}_2\} \le \tilde{\gamma}_1 + \tilde{\gamma}_2$. $\square$

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

As a consequence of the previous results we have the following corollary, which will be useful in Section 2.3:

**Corollary 2.11.** Let $\mu \in M^+(\mathbb{R}^d)$ satisfy (A1) for some $C_\eta \in (0, \infty)$, then for all $\rho \in P(\mathbb{R}^d)$ and $j \in M(G)$ there holds

$$ \frac{1}{2} \iint_G (2 \wedge |x - y|) \eta(x, y) \, d|j|(x, y) \le \sqrt{2C_\eta A(\mu; \rho, j)}. \quad (2.11) $$

**Proof.** Let us consider the case $A(\mu; \rho, j) < \infty$, otherwise the result is trivial. From Lemma 2.6 we have $d\boldsymbol{j}(x, y) = v(x, y)_+ d\gamma_1(x, y) - v(x, y)_- d\gamma_2(x, y)$, with $d\gamma_1(x, y) = d\rho(x)\mu(y)$ and $d\gamma_2(x, y) = d\mu(x)d\rho(y)$. Applying Lemma 2.10 for $\Phi(x, y) = 2 \wedge |x - y|$ and noticing $\Phi(x, y) \le |x - y| \le |x - y| \vee |x - y|^2$, we arrive at the bound

$$ \begin{aligned} & \left(\frac{1}{2} \iint_G (2 \wedge |x - y|) \eta(x, y) \, d\boldsymbol{j}\right)^2 \\ &\le A(\mu; \rho, j) \iint_G (2 \wedge |x - y|)^2 \eta(x, y) (d\gamma_1 + d\gamma_2) \\ &\le A(\mu; \rho, j) 2 \iint_G (|x - y|^2 \vee |x - y|^4) \eta(x, y) \, d\mu(y) \, d\rho(x) \\ &\le A(\mu; \rho, j) 2C_\eta, \end{aligned} $$

where the last estimate follows from (A1) and the integral is finite since $\rho \in P(\mathbb{R}^d)$. $\square$

<figure><img src="image_6.png" /></figure>
**Lemma 2.12.** (Convexity of the action) Let $\mu^i \in \mathcal{M}^+(\mathbb{R}^d)$, $\rho^i \in \mathcal{P}(\mathbb{R}^d)$ and $\mathbf{j}^i \in \mathcal{M}(G)$ for $i = 0, 1$. For $\tau \in (0, 1)$ such that $\mu^\tau = (1 - \tau)\mu^0 + \tau\mu^1$, $\rho^\tau = (1 - \tau)\rho^0 + \tau\rho^1$ and $\mathbf{j}^\tau = (1 - \tau)\mathbf{j}^0 + \tau\mathbf{j}^1$, it holds

$$ \mathcal{A}(\mu^\tau; \rho^\tau, \mathbf{j}^\tau) \leq (1 - \tau)\mathcal{A}(\mu^0; \rho^0, \mathbf{j}^0) + \tau\mathcal{A}(\mu^1; \rho^1, \mathbf{j}^1). $$

**Proof.** Let us consider a measure $\lambda \in \mathcal{M}(G)$ such that $d\gamma_j^i = \tilde{\gamma}_j^i d\lambda$ and $d\mathbf{j}^i = \tilde{\mathbf{j}}^i d\lambda$ for $i = 0, 1$ and $j = 1, 2$. Then, the convex combinations are such that $d\gamma_j^\tau = \tilde{\gamma}_j^\tau d\lambda$ and $d\mathbf{j}^\tau = \tilde{\mathbf{j}}^\tau d\lambda$, where

$$ \tilde{\gamma}_j^\tau = (1 - \tau)\tilde{\gamma}_j^0 + \tau\tilde{\gamma}_j^1, \quad \text{for } j = 1, 2, \\ \text{and} \\ \tilde{\mathbf{j}}^\tau = (1 - \tau)\tilde{\mathbf{j}}^0 + \tau\tilde{\mathbf{j}}^1. $$

Using the convexity of the function $\alpha$ we get the result, that is,

$$ \begin{aligned} \mathcal{A}(\mu^\tau; \rho^\tau, \mathbf{j}^\tau) &= \frac{1}{2} \iint_G \left( \alpha(\tilde{\mathbf{j}}^\tau, \tilde{\gamma}_1^\tau) + \alpha(-\tilde{\mathbf{j}}^\tau, \tilde{\gamma}_2^\tau) \right) \eta(x, y) \, d\lambda(x, y) \\ &\leq \frac{1 - \tau}{2} \iint_G \left( \alpha(\tilde{\mathbf{j}}^0, \tilde{\gamma}_1^0) + \alpha(-\tilde{\mathbf{j}}^0, \tilde{\gamma}_2^0) \right) \eta(x, y) \, d\lambda(x, y) \\ &\quad + \frac{\tau}{2} \iint_G \left( \alpha(\tilde{\mathbf{j}}^1, \tilde{\gamma}_1^1) + \alpha(-\tilde{\mathbf{j}}^1, \tilde{\gamma}_2^1) \right) \eta(x, y) \, d\lambda(x, y) \\ &= (1 - \tau) \mathcal{A}(\mu^0; \rho^0, \mathbf{j}^0) + \tau \mathcal{A}(\mu^1; \rho^1, \mathbf{j}^1). \end{aligned} $$

□

## 2.3. Nonlocal Continuity Equation

In view of the considerations made in Section 2.2, we now deal with the nonlocal continuity equation

$$ \partial_t \rho_t + \bar{\nabla} \cdot \mathbf{j}_t = 0 \quad \text{on } (0, T) \times \mathbb{R}^d, \qquad (2.12) $$

where $(\rho_t)_{t \in [0, T]}$ and $(\mathbf{j}_t)_{t \in [0, T]}$ are unknown Borel families of measures in $\mathcal{P}(\mathbb{R}^d)$ and $\mathcal{M}(G)$, respectively. Equation (2.12) is understood in the weak form: $\forall \varphi \in C_c^\infty((0, T) \times \mathbb{R}^d)$,

$$ \int_0^T \int_{\mathbb{R}^d} \partial_t \varphi_t(x) \, d\rho_t(x) \, dt + \frac{1}{2} \int_0^T \iint_G \bar{\nabla} \varphi_t(x, y) \eta(x, y) \, d\mathbf{j}_t(x, y) \, dt = (2.13) $$

Since $|\bar{\nabla}\varphi(x, y)| \leq ||\varphi||_{C^1}(2 \wedge |x - y|)$, the weak formulation is well-defined under the integrability condition

$$ \int_0^T \iint_G (2 \wedge |x - y|) \eta(x, y) \, d\mathbf{j}_t(x, y) \, dt < \infty. \qquad (2.14) $$

**Remark 2.13.** The integrability condition (2.14) is automatically satisfied by a pair $(\rho_t, \mathbf{j}_t)_{t \in [0, T]}$ such that $\int_0^T \mathcal{A}(\mu; \rho_t, \mathbf{j}_t) \, dt < \infty$, due to Corollary 2.11.
Hence we arrive at the following definition of weak solution of the nonlocal
continuity equation:

**Definition 2.14.** (Nonlocal continuity equation in flux form) A pair $(\rho, j) : [0, T] \to \mathcal{P}(\mathbb{R}^d) \times \mathcal{M}(G)$ is called a weak solution to the nonlocal continuity equation (2.12) provided that

(i) $(\rho_t)_{t \in [0, T]}$ is weakly continuous curve in $\mathcal{P}(\mathbb{R}^d)$;

(ii) $(j_t)_{t \in [0, T]}$ is a Borel-measurable curve in $\mathcal{M}(G)$;

(iii) the pair $(\rho, j)$ satisfies (2.13).

We denote the set of all weak solutions on the time interval $[0, T]$ by $CE_T$. For $\rho^0, \rho^1 \in \mathcal{P}(\mathbb{R}^d)$, a pair $(\rho, j) \in \mathrm{CE}(\rho^0, \rho^1)$ if $(\rho, j) \in \mathrm{CE} := \mathrm{CE}_1$ and in addition $\rho(0) = \rho^0$ and $\rho(1) = \rho^1$.

The following lemma shows that any weak solution satisfying (2.13), which additionally satisfies the integrability condition (2.14) has a weakly continuous representative and hence is a weak solution in the sense of Definition 2.14. This observation justifies the terminology of curve in the space of probability measures; see [2, Lemma 8.1.2] and [23, Lemma 3.1].

**Lemma 2.15.** Let $(\rho_t)_{t \in [0, T]}$ and $(j_t)_{t \in [0, T]}$ be Borel families of measures in $\mathcal{P}(\mathbb{R}^d)$ and $\mathcal{M}(G)$ satisfying (2.13) and (2.14). Then there exists a weakly continuous curve $(\bar{\rho}_t)_{t \in [0, T]} \subset \mathcal{P}(\mathbb{R}^d)$ such that $\bar{\rho}_t = \rho_t$ for a.e. $t \in [0, T]$. Moreover, for any $\varphi \in C_c^\infty([0, T] \times \mathbb{R}^d)$ and all $0 \le t_0 \le t_1 \le T$ it holds that

$$
\begin{align}
& \int_{\mathbb{R}^d} \varphi_{t_1}(x) \, d\bar{\rho}_{t_1}(x) - \int_{\mathbb{R}^d} \varphi_{t_0}(x) \, d\bar{\rho}_{t_0}(x) \nonumber \\
& = \int_{t_0}^{t_1} \int_{\mathbb{R}^d} \partial_t \varphi_t(x) \, d\rho_t(x) \, dt \tag{2.15} \\
& \quad + \frac{1}{2} \int_{t_0}^{t_1} \iint_G \bar{\nabla} \varphi_t(x, y) \eta(x, y) \, d\mathbf{j}_t(x, y) \, dt. \nonumber
\end{align}
$$

We now prove propagation of second-order moments.

**Lemma 2.16.** (Uniformly bounded second moments) Let $(\mu^n)_n \subset \mathcal{M}^+(\mathbb{R}^d)$ such that (A1) holds uniformly in $n$. Let $(\rho_0^n)_n \subset \mathcal{P}_2(\mathbb{R}^d)$ be such that $\sup_{n \in \mathbb{N}} M_2(\rho_0^n) < \infty$ and $(\rho^n, j^n)_n \subset \mathrm{CE}_T$ be such that $\sup_{n \in \mathbb{N}} \int_0^T \mathcal{A}(\mu^n; \rho_t^n, j_t^n) dt < \infty$. Then $\sup_{t \in [0, T]} \sup_{n \in \mathbb{N}} M_2(\rho_t^n) < \infty$.

**Proof.** We proceed by considering the time derivative of the second-order moment of $\rho_t^n$ for all $t \in [0, T]$ and $n \in \mathbb{N}$. Since $x \mapsto |x|^2$ is not an admissible test function in (2.13), we introduce a smooth cut-off function $\varphi_R$ satisfying $\varphi_R(x) = 1$ for $x \in B_R$, $\varphi_R(x) = 0$ for $x \in \mathbb{R}^d \setminus B_{2R}$ and $|\nabla \varphi_R| \le \frac{2}{R}$. Then, we can use the definition of solution with the function $\psi_R(x) = \varphi_R(x)^2(|x|^2 + 1)$ and apply Lemma 2.10 with $\Phi = \bar{\nabla}\psi_R$ to obtain, for all $t \in [0, T]$ and $n \in \mathbb{N}$,

$$
\frac{d}{dt} \int_{\mathbb{R}^d} \psi_R(x) \, d\rho_t^n(x)
$$
$$ = \frac{1}{2} \iint_G \bar{\nabla} \psi_R(x, y) \eta(x, y) \, d\mathbf{j}_t^n(x, y) \\ \leq \sqrt{\mathcal{A}(\mu^n; \rho_t^n, \mathbf{j}_t^n)} \left( \iint_G |\bar{\nabla} \psi_R(x, y)|^2 \eta(x, y) (d\gamma_1^n + d\gamma_2^n) \right)^{\frac{1}{2}}. $$

For $R \ge 1$, we estimate, for all $(x, y) \in G$,

$$ |\bar{\nabla} \psi_R(x, y)|^2 \le 2|\varphi_R(y)^2 - \varphi_R(x)^2|^2 + 2|\varphi_R(y)^2|y^2 - \varphi_R(x)^2|x|^2, \quad (2.16) $$

and observe that

$$ |\bar{\nabla} \varphi_R^2(x, y)| = |\bar{\nabla} \varphi_R(x, y) (\varphi_R(x) + \varphi_R(y))| \le \frac{4}{R} |x - y|. $$

Hence the first term in (2.16) is bounded by $32|x - y|^2$, since $R \ge 1$. For the second term in (2.16), we abbreviate by setting $r = \varphi_R(x)|x|$ and $s = \varphi_R(y)|y|$ and compute the bound

$$ \begin{aligned} |s^2 - r^2|^2 &= |s - r|^2 |s + r|^2 \le 2|s - r|^4 + 8|r|^2 |s - r|^2 \\ &\le 8 (|r|^2 + 1) (|s - r|^2 \lor |s - r|^4). \end{aligned} $$

It is easy to check that $x \mapsto \varphi_R(x)|x|$ is globally Lipschitz and we can conclude that, for some numerical constant $C > 0$, for all $(x, y) \in G$ we have

$$ \begin{aligned} |\bar{\nabla} \psi_R(x, y)|^2 &\le 32|x - y|^2 + C|x|^2 (|x - y|^2 \lor |x - y|^4) \\ &\le C (|x|^2 + 1) (|x - y|^2 \lor |x - y|^4). \end{aligned} $$

Thus, by sending $R \to \infty$ and using (A1), it follows that

$$ \frac{d}{dt} \int_{\mathbb{R}^d} (|x|^2 + 1) \, d\rho_t^n(x) \le \sqrt{\mathcal{A}(\mu^n; \rho_t^n, \mathbf{j}_t^n)} \left( 2CC_\eta \int_{\mathbb{R}^d} (|x|^2 + 1) \, d\rho_t^n(x) \right)^{\frac{1}{2}} $$

By integrating the above differential inequality, we arrive at the bound

$$ \int_{\mathbb{R}^d} |x|^2 \, d\rho_t^n(x) \le 2 \int_{\mathbb{R}^d} (|x|^2 + 1) \, d\rho_0^n(x) + 2CC_\eta T \int_0^T \mathcal{A}(\mu^n; \rho_t^n, \mathbf{j}_t^n) \, dt, $$

whence we conclude by taking the suprema in $n \in \mathbb{N}$ and $t \in [0, T]$. $\square$

Now we are ready to show compactness for the solutions to (2.12).

**Proposition 2.17.** (Compactness of solutions to the nonlocal continuity equation) Let $(\mu^n)_n \subset \mathcal{M}^+(\mathbb{R}^d)$ and suppose that $(\mu^n)_n$ narrowly converges to $\mu$. Moreover, suppose that the base measures $\mu^n$ and $\mu$ satisfy (A1) and (A2) uniformly in $n$. Let $(\rho^n, \mathbf{j}^n) \in \mathrm{CE}_T$ for each $n \in \mathbb{N}$ be such that $(\rho_0^n)_n$ satisfies $\sup_{n \in \mathbb{N}} M_2(\rho_0^n) < \infty$ and

$$ \sup_{n \in \mathbb{N}} \int_0^T \mathcal{A}(\mu^n; \rho_t^n, \mathbf{j}_t^n) \, dt < \infty. \quad (2.17) $$
Then, there exists $(\rho, \mathbf{j}) \in \mathrm{CE}_T$ such that, up to a subsequence, as $n \to \infty$ it holds

$$ \begin{aligned} & \rho_t^n \to \rho_t \quad \text{for all } t \in [0, T], \\ & \mathbf{j}^n \to \mathbf{j} \quad \text{in } \mathcal{M}_{\mathrm{loc}}(G \times [0, T]), \end{aligned} $$

with $\rho_t \in \mathcal{P}_2(\mathbb{R}^d)$ for any $t \in [0, T]$. Moreover, the action is lower semicontinuous along the above subsequences $(\mu^n)_n$, $(\rho^n)_n$ and $(\mathbf{j}^n)_n$, i.e.,

$$ \liminf_{n \to \infty} \int_0^T \mathcal{A}(\mu^n; \rho_t^n, \mathbf{j}_t^n) \, dt \ge \int_0^T \mathcal{A}(\mu; \rho_t, \mathbf{j}_t) \, dt. $$

**Proof.** We argue similarly to [22, Lemma 4.5], [23, Proposition 3.4]. For each $n \in \mathbb{N}$ we define $\mathbf{j}^n \in \mathcal{M}(G \times [0, T])$ as $d\mathbf{j}^n(x, y, t) = d\mathbf{j}_t^n(x, y) dt$. In view of Lemma 2.16 there exists $C_2 > 0$ such that $\sup_{t \in [0, T]} \sup_{n \in \mathbb{N}} M_2(\rho_t^n) \le C_2 < +\infty$.

For any compact sets $K \subset G$ and $I \subseteq [0, T]$, we apply the bound (2.11) of Corollary 2.11 and the Cauchy–Schwarz inequality to get

$$ \begin{aligned} \sup_{n \in \mathbb{N}} |\mathbf{j}^n|(K \times I) &\le \sup_{n \in \mathbb{N}} \int_I \iint_K \frac{(2 \wedge |x-y|) \eta(x, y)}{(2 \wedge |x-y|) \eta(x, y)} d|\mathbf{j}_t^n|(x, y) dt \\ &\le \frac{2\sqrt{|I|}\sqrt{2C_\eta}}{\inf_{(x,y) \in K} (2 \wedge |x-y|)\eta(x,y)} \left( \sup_{n \in \mathbb{N}} \int_0^T \mathcal{A}(\mu^n; \rho_t^n, \mathbf{j}_t^n) dt \right)^{\frac{1}{2}}. \end{aligned} \quad (2.18) $$

Thanks to Assumption (W), we have that $\inf_{(x,y) \in K} (2 \wedge |x-y|)\eta(x, y) > 0$ for any compact $K \subset G$. Hence, by (2.17), $(\mathbf{j}^n)_n$ has total variation uniformly bounded in $n$ on every compact set of $G \times [0, T]$, which implies, up to a subsequence, $\mathbf{j}^n \to \mathbf{j}$ as $n \to \infty$ in $\mathcal{M}_{\mathrm{loc}}(G \times [0, T])$. Because of the disintegration theorem, there exists a Borel family $(\mathbf{j}_t)_{t \in [0, T]}$ such that, for all compact sets $I \subseteq [0, T]$ and $K \subset G$, there holds that $\mathbf{j}(K \times I) = \int_I \mathbf{j}_t(K) dt$. Thanks to the bound (2.18), the family $\{\mathbf{j}_t\}_{t \in [0, T]}$ still satisfies (2.14).

Now, as we need to pass to the limit in (2.13), we consider a function $\xi \in C_c^\infty(\mathbb{R}^d)$ and an interval $[t_0, t_1] \subseteq [0, T]$. The function $\chi_{[t_0, t_1]}(t)\bar{\nabla}\xi(x, y)$ has no compact support in $[t_0, t_1] \times G$, so we proceed by a truncation argument. Let $\varepsilon > 0$ and let us set $I^\varepsilon = [t_0 + \varepsilon, t_1 - \varepsilon]$, $N_\varepsilon = \bar{B}_{\varepsilon^{-1}} \times \bar{B}_{\varepsilon^{-1}}$, where $B_{\varepsilon^{-1}} = \{x \in \mathbb{R}^d : |x| < \varepsilon^{-1}\}$, and $G_\varepsilon = \{(x, y) \in G : \varepsilon \le |x-y|\}$. Hence we can find $\varphi_\varepsilon \in C_c^\infty([t_0, t_1] \times G; [0, 1])$ satisfying

$$ \{\varphi_\varepsilon = 1\} \supseteq I_\varepsilon \times (G_\varepsilon \cap N_\varepsilon), \quad (2.19) $$

so that $\varphi_\varepsilon \to \chi_{[t_0, t_1]} \chi_G$ as $\varepsilon \to 0$ and $\varphi_\varepsilon \chi_{[t_0, t_1]} \bar{\nabla}\xi$ has compact support in $[t_0, t_1] \times G$. Then, we get thanks to Assumption (W), that

$$ \begin{aligned} & \lim_{n \to \infty} \int_{t_0}^{t_1} \iint_G \varphi_\varepsilon(t, x, y) \bar{\nabla} \xi(x, y) \eta(x, y) d\mathbf{j}_t^n(x, y) dt \\ = & \int_{t_0}^{t_1} \iint_G \varphi_\varepsilon(t, x, y) \bar{\nabla} \xi(x, y) \eta(x, y) d\mathbf{j}_t(x, y) dt. \end{aligned} \quad (2.20) $$
Now, it remains to show that

$$ \lim_{\epsilon \to 0} \sup_{n \in \mathbb{N}} \left| \int_{t_0}^{t_1} \iint_G (1 - \varphi_{\epsilon}(t, x, y)) \bar{\nabla} \xi(x, y) \eta(x, y) \, d\mathbf{j}_t^n(x, y) \, dt \right| = 0. \quad (2.21) $$

We need to estimate terms for which $\varphi_{\epsilon}(t, x) < 1$. First, setting $I_{\epsilon}^{c} = [t_0, t_1] \setminus I_{\epsilon}$, we note that

$$ [t_0, t_1] \times G \setminus \{\varphi_{\epsilon} = 1\} \subseteq (I_{\epsilon}^{c} \times G) \cup (I_{\epsilon} \times (G \setminus (G_{\epsilon} \cap N_{\epsilon}))) =: M_{\epsilon}, $$

whence, by Lemma 2.10,

$$ \begin{aligned} & \left| \int_{t_0}^{t_1} \iint_G (1 - \varphi_{\epsilon}(t, x, y)) \bar{\nabla} \xi(x, y) \eta(x, y) \, d\mathbf{j}_t^n(x, y) \, dt \right| \\ & \le \|\xi\|_{C^1} \int_{t_0}^{t_1} \iint_G (1 - \varphi_{\epsilon}(t, x, y)) (2 \wedge |x - y|) \eta(x, y) \, d|\mathbf{j}_t^n|(x, y) \, dt \\ & \le 2 \|\xi\|_{C^1} \left( \int_0^T \mathcal{A}(\mu^n; \rho_t^n, \mathbf{j}_t^n) \, dt \right)^{\frac{1}{2}} \\ & \quad \times \left( \iiint_{M_{\epsilon}} (4 \wedge |x - y|^2) \eta(x, y) \, d(\gamma_{1,t}^n + \gamma_{2,t}^n) \, dt \right)^{\frac{1}{2}}. \end{aligned} $$

Since $4 \wedge |x - y|^2 \le |x - y|^2 \vee |x - y|^4$ we have, by Assumption (A1), the bound

$$ \int_{I_{\epsilon}^{c}} \iint_{G} (4 \wedge |x - y|^2) \eta(x, y) \, d(\gamma_{t}^{1,n} + \gamma_{t}^{2,n}) \, dt \le 2|I_{\epsilon}^{c}|C_{\eta} = 4C_{\eta}\epsilon. $$

Likewise, using the symmetry, we arrive at

$$ \begin{aligned} & \int_{I_{\epsilon}} \iint_{G_{\epsilon}^{c}} (4 \wedge |x - y|^2) \eta(x, y) \, d(\gamma_{1,t}^{n} + \gamma_{2,t}^{n}) \, dt \\ & = 2 \int_{0}^{T} \iint_{G_{\epsilon}^{c}} (4 \wedge |x - y|^2) \eta(x, y) \, d\mu^{n}(y) \, d\rho_{t}^{n}(x) \, dt, \end{aligned} $$

which vanishes as $\epsilon \to 0$ in view of Assumption (A2). Finally, the last term is estimated again using (A1):

$$ \begin{aligned} & \int_{I_{\epsilon}} \iint_{G \setminus N_{\epsilon}} (4 \wedge |x - y|^2) \eta(x, y) \, d\gamma_t^{1,n} \, dt \\ & \le \int_0^T \iint_{\bar{B}_{\epsilon-1}^c} \int_{\mathbb{R}^d} (4 \wedge |x - y|^2) \eta(x, y) \, d\mu^n(y) \, d\rho_t^n(x) \, dt \\ & \le T C_{\eta} \sup_{t \in [0, T]} \rho_t^n(\bar{B}_{\epsilon-1}^c) \to 0 \quad \text{as } \epsilon \to 0, \end{aligned} $$

since $M_2(\rho_t^n) \le C_2$ for any $n \in \mathbb{N}$ and $t \in [0, T]$ by Lemma 2.16.
Combining (2.20) and (2.21), we get

$$ 
\begin{aligned}
& \lim_{n \to \infty} \int_{t_0}^{t_1} \iint_G \bar{\nabla} \xi(x, y) \eta(x, y) \, d\mathbf{j}_t^n(x, y) \, dt \\
& = \int_{t_0}^{t_1} \iint_G \bar{\nabla} \xi(x, y) \eta(x, y) \, d\mathbf{j}_t(x, y) \, dt.
\end{aligned}
$$

By means of the last convergence, the tightness of $(\rho_0^n)_n$, and (2.15) with $\varphi(t, x) = \xi(x)$, $t_0 = 0$ and $t_1 = T$, we obtain that $(\rho_t^n)_n$ locally narrowly converges to some finite non-negative measure $\rho_t \in M^+(\mathbb{R}^d)$ for any $t \in [0, T]$. In particular, for any $\xi \in C_c^\infty(\mathbb{R}^d)$ and any $t \in [0, T]$, we have

$$ \int_{\mathbb{R}^d} \xi(x) \, d\rho_t(x) = \int_{\mathbb{R}^d} \xi(x) \, d\rho_0(x) + \frac{1}{2} \int_0^t \iint_G \bar{\nabla} \xi(x, y) \eta(x, y) \, d\mathbf{j}_s(x, y) \, ds. $$

Now, for $R > 0$, let us consider a function $\xi_R \in C_c^\infty(\mathbb{R}^d)$ such that $0 \le \xi \le 1$, $\xi = 1$ on $B_R$, and $\|\xi\|_{C^1} \le 1$. Because of the integrability condition (2.14), satisfied thanks to Corollary 2.11, we have

$$ 
\begin{aligned}
& \left| \int_0^t \frac{1}{2} \iint_G \bar{\nabla} \xi_R(x, y) \eta(x, y) \, d\mathbf{j}_s(x, y) \, ds \right| \\
& \le \frac{1}{2} \int_0^t \iint_{G \setminus (B_R \times B_R)} (2 \wedge |x - y|) \eta(x, y) \, d|\mathbf{j}_s| \, ds \xrightarrow{R \to \infty} 0.
\end{aligned}
$$

Hence the measure $\rho_t$ is actually a probability measure on $\mathbb{R}^d$ for all $t \in [0, T]$. Moreover Lemma 2.16 ensures that the convergence is global and not only local. As a direct consequence of the previous considerations, $(\rho, \mathbf{j}) \in CE_T$ and the lower semicontinuity follows from Lemma 2.9. □

## 2.4. Nonlocal Upwind Transportation Quasi-Metric

Here, we give a rigorous definition of the nonlocal transportation quasi-metric we introduced in (1.8). Let us recall that $\eta: \{(x, y) \in \mathbb{R}^d \times \mathbb{R}^d : x \ne y\} \to [0, \infty)$ is the weight function satisfying (W).

**Definition 2.18.** (Nonlocal upwind transportation cost) For $\mu \in M^+(\mathbb{R}^d)$ satisfying Assumptions (A1) and (A2), and $\rho_0, \rho_1 \in P_2(\mathbb{R}^d)$, the nonlocal upwind transportation cost between $\rho_0$ and $\rho_1$ is defined by

$$ T_{\mu}(\rho_0, \rho_1)^2 = \inf \left\{ \int_0^1 A(\mu; \rho_t, \mathbf{j}_t) \, dt : (\rho, \mathbf{j}) \in CE(\rho_0, \rho_1) \right\}. \quad (2.22) $$

If $\mu$ is clear from the context, the notation $T$ is used in place of $T_\mu$.

Note that Proposition 2.17 ensures the existence of minimizers to (2.22), when $T_\mu < \infty$, which holds when there exists a path of finite action. On the other hand, if this is not the case, the nonlocal upwind transportation cost is infinite. For example, consider the graph with vertices set by $\mu$ and $\eta$ which is disconnected, meaning that
there are $x, y \in \text{supp } \mu$ such that there is no sequence $(x_0 = x, x_1, \dots, x_{n-1}, x_n = y)_n$ with $\eta(x_i, x_{i+1}) > 0$ for all $i = 0, \dots, n-1$; in this case, $T_\mu(\delta_x, \delta_y) = \infty$ since the set of solutions to the continuity equation $\text{CE}(\delta_x, \delta_y)$ is empty.

Due to the one-homogeneity of the action density function $\alpha$ in (2.5), we have the following reparametrization result, which is similar to [22, Theorem 5.4]:

**Lemma 2.19.** (Reparametrization) For any $\mu \in \mathcal{M}^+(\mathbb{R}^d)$ satisfying Assumptions (A1) and (A2), and any $\rho_0, \rho_T \in \mathcal{P}_2(\mathbb{R}^d)$, it holds that

$$ T_{\mu}(\rho_0, \rho_T) = \inf \left\{ \int_{0}^{T} \sqrt{A(\mu; \rho_t, j_t)} \, dt : (\rho, j) \in \text{CE}_T(\rho_0, \rho_T) \right\}. $$

Now, as consequence of the above reparametrization and Jensen's inequality, we have the following result, which implies that the infimum is in fact a minimum; see [23, Proposition 4.3].

**Proposition 2.20.** For any $\mu \in \mathcal{M}^+(\mathbb{R}^d)$ satisfying Assumptions (A1) and (A2), and any $\rho_0, \rho_1 \in \mathcal{P}_2(\mathbb{R}^d)$ such that $T_\mu(\rho_0, \rho_1) < \infty$, the infimum in (2.22) is attained by a curve $(\rho, j) \in \text{CE}(\rho_0, \rho_1)$ so that $A(\rho_t, j_t) = T_\mu(\rho_0, \rho_1)^2$ for a.e. $t \in [0, 1]$. Such curve is a constant-speed geodesic for $T_\mu$, i.e.,

$$ T_{\mu}(\rho_s, \rho_t) = |t - s| T_{\mu}(\rho_0, \rho_1), \quad \text{for all } s, t \in [0, 1]. $$

The next proposition establishes a link between $T_\mu$ and the $W_1$-distance.

**Proposition 2.21.** (Comparison with $W_1$) Let $\mu \in \mathcal{M}^+(\mathbb{R}^d)$ satisfy (A1) for some $C_\eta > 0$ (depending only on $\mu$ and $\eta$). Then for any $\rho^0, \rho^1 \in \mathcal{P}_2(\mathbb{R}^d)$ it holds

$$ W_1(\rho^0, \rho^1) \le \sqrt{2C_\eta} \sqrt{T(\rho^0, \rho^1)}. $$

**Proof.** By a standard regularization argument and the truncation procedure as in the proof of Lemma 2.16, we can actually consider any 1-Lipschitz function $\psi$ as a test function in the weak formulation (2.13) for some $(\rho, j) \in \text{CE}(\rho^0, \rho^1)$. Then we can estimate, by Lemma 2.10 and Assumption (A1),

$$
\begin{aligned}
& \left| \int_{\mathbb{R}^d} \psi \, d\rho^1 - \int_{\mathbb{R}^d} \psi \, d\rho^0 \right| \\
&= \left| \frac{1}{2} \int_0^1 \iint_G \bar{\nabla} \psi \, \eta \, d\mathbf{j}_t \, dt \right| \le \frac{1}{2} \int_0^1 \iint_G |x-y| \eta(x, y) \, d|\mathbf{j}_t|(x, y) \, dt \\
&\le \left( \int_0^1 A(\rho_t, \mathbf{j}_t) \, dt \right)^{\frac{1}{2}} \left( \int_0^1 \iint_G |x-y|^2 \eta(x, y) (d\gamma_1 + d\gamma_2) \right)^{\frac{1}{2}} \\
&\le \left( \int_0^1 A(\rho_t, \mathbf{j}_t) \, dt \right)^{\frac{1}{2}} \\
&\quad \times \left( 2 \int_0^1 \iint_G (|x-y|^2 \vee |x-y|^4) \eta(x, y) \, d\mu(y) \, d\rho_t(x) \right)^{\frac{1}{2}}
\end{aligned}
$$
$$ \leq \sqrt{2C_{\eta}} \left( \int_{0}^{1} A(\rho_{t}, j_{t}) dt \right)^{\frac{1}{2}}. $$

Taking the supremum over all 1-Lipschitz functions and the infimum in the couplings $(\rho, j) \in \text{CE}(\rho^0, \rho^1)$ gives the result. □

The results above show that $\mathcal{T}_{\mu}$ is an extended (meaning that it can take value $\infty$) quasi-metric on the set of probability measures which induces a topology stronger than the $W_1$-topology.

**Theorem 2.22.** Let $\mu \in \mathcal{M}^{+}(\mathbb{R}^d)$ satisfy Assumptions (A1) and (A2). The nonlocal upwind transportation cost $\mathcal{T}_{\mu}$ defines an extended quasi-metric on $\mathcal{P}_2(\mathbb{R}^d)$. The map $(\rho_0, \rho_1) \mapsto \mathcal{T}_{\mu}(\rho_0, \rho_1)$ is lower semicontinuous with respect to the narrow convergence. The topology induced by $\mathcal{T}_{\mu}$ is stronger than the $W_1$-topology and the narrow topology. In particular, bounded sets are narrowly relatively compact in $(\mathcal{P}_2(\mathbb{R}^d), \mathcal{T}_{\mu})$.

**Proof.** If $\mathcal{T}_{\mu}(\rho_0, \rho_1) = 0$, then $\mathcal{A}(\mu; \rho_t, j_t) = 0$ for a.e. $t \in [0, 1]$. Hence $j_t \equiv 0$ $\gamma_t$-a.e., which implies that $\rho_0 \equiv \rho_1$ by the nonlocal continuity equation (2.15). The triangle inequality is a consequence of Lemma 2.19 and the fact that solutions to the nonlocal continuity equation can be concatenated. The lower semicontinuity and compactness properties of $\mathcal{T}_{\mu}$ are inherited from the action functional $\mathcal{A}$ via Proposition 2.17. In view of the comparison with $W_1$ from Proposition 2.21, we have that the topology induced by $\mathcal{T}_{\mu}$ is stronger than that induced by $W_1$ and the narrow topology. □

The next lemma provides a quantitative illustration of asymmetry of $\mathcal{T}$.

**Lemma 2.23.** (Two-point space) Let us consider the two-point graph $\Omega := \{0, 1\}$, with $\eta(0, 1) = \eta(1, 0) = \alpha > 0$, $\mu(0) = p > 0$ and $\mu(1) = q > 0$. Let $\rho, v \in \mathcal{P}_2(\Omega)$ and let $\rho_0, \rho_1, v_0, v_1 \in [0, 1]$ be such that $\rho = \rho_0\delta_0 + \rho_1\delta_1$ and $v = v_0\delta_0 + v_1\delta_1$. There holds

$$ \mathcal{T}(\rho, v) = \begin{cases} \frac{2}{\sqrt{\alpha p}} (\sqrt{\rho_1} - \sqrt{v_1}) & \text{if } \rho_0 < v_0, \\ \frac{2}{\sqrt{\alpha q}} (\sqrt{\rho_0} - \sqrt{v_0}) & \text{if } v_0 < \rho_0. \end{cases} \quad (2.23) $$

**Proof.** Let us fix $\lambda = \delta_{(0,1)} + \delta_{(1,0)}$ and notice that $\rho_0 + \rho_1 = 1$ and $v_0 + v_1 = 1$ as $\rho, v$ are probability measures. Since $\Omega = \{0, 1\}$, note that for any curve $t \in [0, 1] \mapsto \rho_t \in \mathcal{P}_2(\Omega)$ there exists a function $g: t \in [0, 1] \mapsto g_t \in [0, 1]$ accounting for the mass displacement. Thus, we notice that $(\rho, j) \in \text{CE}(\rho, v)$ if

$$ \rho_t = g_t \delta_0 + (1 - g_t) \delta_1, \quad j_t \text{ such that } j_t(0, 1) = -\frac{\dot{g}_t}{\alpha} \text{ and } j_t(1, 0) = \frac{\dot{g}_t}{\alpha}, \\ \text{for all } t \in [0, 1]. $$

Hence, using that $j_t$ is antisymmetric yields

$$ \mathcal{T}(\rho_0, \rho_1)^2 = \inf \left\{ \int_0^1 A(\rho_t, j_t) dt : (\rho, j) \in \text{CE}(\rho_0, \rho_1) \right\} $$
<figure><img src="image_7.png" /><figcaption>**Fig. 3.** In the context of Lemma 2.23, in the above figures we parametrize the quasi-distance $\mathcal{T}(\rho, \nu)$ by $\rho_0 \in [0, 1]$ and $\nu_0 \in [0, 1]$ for $\mu(0) = 0.1$ (left) and $\mu(1) = 0.5$ (right). Colors represent different values of $\mathcal{T}(\rho, \nu)$ with respect to the initial values $\rho_0$ and $\nu_0$. In the left figure, by swapping the values of $\rho_0$ and $\nu_0$ on the axes, we can see that $\mathcal{T}$ is non-symmetric</figcaption></figure>

$$ = \inf_g \left\{ \int_0^1 \frac{|(\dot{g}_t)|^2}{\alpha g_t q} + \frac{|(\dot{g}_t)_+|^2}{\alpha(1-g_t)p} \, dt : g_0 = \rho_0 \text{ and } g_1 = \nu_0 \right\}. $$

Now, let us assume without loss of generality that $\rho_0 < \nu_0$. Obviously, in this configuration we can restrict the above infimum among non-decreasing $g$, as it gives a lower action. Therefore, by applying Jensen's inequality, we have

$$ \begin{aligned} \mathcal{T}(\rho, \nu)^2 &= \inf_{g \nearrow} \frac{1}{\alpha p} \int_0^1 \frac{|\dot{g}_t|^2}{(1-g_t)} \, dt = \inf_{g \nearrow} \frac{1}{\alpha p} \int_0^1 \left| -2 \frac{d}{dt} \left( \sqrt{1-g_t} \right) \right|^2 \, dt \\ &\geq \inf_{g \nearrow} \frac{4}{\alpha p} \left| \int_0^1 -\frac{d}{dt} \left( \sqrt{1-g_t} \right) \, dt \right|^2 = \frac{4}{\alpha p} \left( \sqrt{1-\rho_0} - \sqrt{1-\nu_0} \right)^2 \\ &= \frac{4}{\alpha p} \left( \sqrt{\rho_1} - \sqrt{\nu_1} \right)^2. \end{aligned} $$

The equality case is obtained by noting that the solution to $-\frac{d}{dt}\sqrt{1-g_t} = \sqrt{\rho_1} - \sqrt{\nu_1}$ for all $t \in [0, 1]$, with consistent boundary values $g_0 = \rho_0$ and $g_1 = \nu_0$, is given by $g_t = 1 - (\sqrt{\rho_1}(1-t) + \sqrt{\nu_1}t)^2$. The case $\nu_0 < \rho_0$ is obtained in a similar manner, which gives formula (2.23). □

**Remark 2.24.** The quasi-metric is in general already non-symmetric on the two-point space, which one can best observe in Fig. 3. In the case $p = \frac{1}{2}$, the swapping $\hat{\rho}_0 = \rho_1$ and $\hat{\rho}_1 = \rho_0$ preserves the quasi-distance $\mathcal{T}(\rho, \nu) = \mathcal{T}(\hat{\rho}, \hat{\nu})$.

We now adapt the standard definition of absolutely continuous curves in metric spaces from [2, Chapter 1] to our setting. Let $\mu \in \mathcal{M}^+(\mathbb{R}^d)$ satisfy Assumptions (A1) and (A2). A curve $[0, T] \ni t \mapsto \rho_t \in \mathcal{P}_2(\mathbb{R}^d)$ is said to be 2-absolutely continuous with respect to $\mathcal{T}_\mu$ if there exists $m \in L^2((0, T))$ such that

$$ \mathcal{T}_\mu(\rho_{t_0}, \rho_{t_1}) \leq \int_{t_0}^{t_1} m(t) \, dt \quad \text{for all } 0 < t_0 \leq t_1 < T. \qquad (2.24) $$
In this case, we write $\rho \in AC([0, T]; (\mathcal{P}_2(\mathbb{R}^d), \mathcal{T}_\mu))$. For any $\rho \in AC([0, T]; (\mathcal{P}_2(\mathbb{R}^d), \mathcal{T}_\mu))$ the quantity

$$|\rho'_t| := \lim_{h \to 0} \frac{\mathcal{T}_\mu(\rho_t, \rho_{t+h})}{|h|} \quad (2.25)$$

is well-defined for a.e. $t \in [0, T]$ and is called the *metric derivative* of $\rho$ at $t$. Moreover, the function $t \to |\rho'(t)|$ belongs to $L^2((0, T))$ and it satisfies $|\rho'(t)| \le m(t)$ for a.e. $t \in [0, T]$, which means $\rho'$ is the minimal integrand satisfying (2.24). The length of a curve $\rho \in AC([0, T]; (\mathcal{P}_2(\mathbb{R}^d), \mathcal{T}_\mu))$ is defined by $L(\rho) := \int_0^T |\rho'(t)| dt$.

**Proposition 2.25.** (Metric velocity) Let $\mu \in \mathcal{M}^+(\mathbb{R}^d)$ satisfy Assumptions (A1) and (A2). A curve $(\rho_t)_{t \in [0,T]} \subset \mathcal{P}_2(\mathbb{R}^d)$ belongs to $AC([0, T]; (\mathcal{P}_2(\mathbb{R}^d), \mathcal{T}_\mu))$ if and only if there exists a family $(\mathbf{j}_t)_{t \in [0,T]}$ such that $(\rho, \mathbf{j}) \in \mathrm{CE}_T$ and

$$\int_0^T \sqrt{\mathcal{A}(\mu; \rho_t, \mathbf{j}_t)} \, dt < \infty.$$

In this case, the metric derivative is bounded as in $|\rho'|^2(t) \le \mathcal{A}(\mu; \rho_t, \mathbf{j}_t)$ for a.e. $t \in [0, T]$. In addition, there exists a unique family $(\mathbf{j}_t)_{t \in [0,T]}$ such that $(\rho, \tilde{\mathbf{j}}) \in \mathrm{CE}_T$ and

$$|\rho'|^2(t) = \mathcal{A}(\mu; \rho_t, \tilde{\mathbf{j}}_t) \quad \text{for a.e. } t \in [0, T]. \qquad (2.26)$$

Hereby, the previous identity holds if and only if $\tilde{\mathbf{j}}_t \in T_\rho \mathcal{P}_2(\mathbb{R}^d)$ for a.e. $t \in [0, T]$, where

$$T_\rho \mathcal{P}_2(\mathbb{R}^d) = \{ \mathbf{j} \in \mathcal{M}_{\gamma_1}^{\text{as}}(G) : \mathcal{A}(\mu; \rho, \mathbf{j}) < \infty, \\ \mathcal{A}(\mu; \rho, \mathbf{j}) \le \mathcal{A}(\mu; \rho, \mathbf{j} + \mathbf{d}) \text{ for all } \mathbf{d} \in \mathcal{M}_{\text{div}}(G) \}, \quad (2.27)$$

with $\mathcal{M}_{\gamma_1}^{\text{as}}(G)$ defined in (2.9), and $\mathcal{M}_{\text{div}}(G)$ the set of nonlocal divergence-free fluxes, that is

$$\mathcal{M}_{\text{div}}(G) = \left\{ \mathbf{d} \in \mathcal{M}(G) : \iint_G \bar{\nabla} \psi \, \eta \, d\mathbf{d} = 0 \text{ for all } \psi \in C_c^\infty(\mathbb{R}^d) \right\}.$$

**Proof.** The first statement on the characterization of absolutely continuous curves as curves of finite action follows from [22, Theorem 5.17], in view of Lemma 2.19 and Propositions 2.17 and 2.20. Let us now show that (2.26) holds if and only if $\tilde{\mathbf{j}}_t$ belongs to $T_\rho \mathcal{P}_2(\mathbb{R}^d)$ for a.e. $t \in [0, 1]$, given by (2.27). Let $t \in [0, 1]$ be so that $\mathbf{j}_t$ verifies $\mathcal{A}(\mu; \rho_t, \mathbf{j}_t) < +\infty$. Due to Corollary 2.8, the element $\tilde{\mathbf{j}}_t$ of minimal action satisfying (2.26) is characterized by $\partial_t \rho_t + \bar{\nabla} \cdot \mathbf{j}_t = 0 = \partial_t \rho_t + \bar{\nabla} \cdot \tilde{\mathbf{j}}_t$, that is,

$$\tilde{\mathbf{j}}_t = \arg \min_{j \in \mathcal{M}_{\gamma_1}^{\text{as}}(G)} \{ \mathcal{A}(\mu; \rho_t, j) : \bar{\nabla} \cdot j = \bar{\nabla} \cdot \tilde{\mathbf{j}}_t \}.$$

Recalling the notation for the Jordan decomposition of a measure from Section 2.2, note that we use that the functional $j \mapsto \mathcal{A}(\mu; \rho, j)$ is strictly convex for $j \in$
$\mathcal{M}(G)$ such that $j^+ \ll \rho \otimes \mu$ and $j^- \ll \mu \otimes \rho$, which is guaranteed above since $\mathcal{A}(\mu; \rho, j) < \infty$ and $j \in \mathcal{M}_{\gamma_1}^{\text{as}}(G)$. Then, we observe the set $\{j \in \mathcal{M}_{\gamma_1}^{\text{as}}(G) : \bar{\nabla} \cdot j = \bar{\nabla} \cdot j_t\}$ is closed with respect to the narrow convergence. In addition, the estimate (2.10) from Lemma 2.10 with $\Phi(x, y) = |x - y| \vee |x - y|^2$ gives

$$ \frac{1}{2} \iint_K \eta(x, y) \, d|j|(x, y) \le \frac{\sqrt{2C_\eta} \sqrt{\mathcal{A}(\mu; \rho_t, j)}}{\inf_K (|x-y| \vee |x-y|^2)} \quad \text{for all compact } K \subset G, $$

showing that the sublevel sets of $j \mapsto \mathcal{A}(\mu; \rho_t, j)$ are locally relatively compact with respect to the narrow convergence, arguing as in the proof of Proposition 2.17. Hence the element $\tilde{j}_t$ is well-defined by applying the direct method of calculus of variations. $\square$

We defined the tangent space $T_\rho P_2(\mathbb{R}^d)$ in (2.27) using the nonlocal fluxes $j$. We note that this is in some way a nonlocal, Lagrangian description of the tangent vectors and that the relationship between this Lagrangian description and the Eulerian description is the nonlocal continuity equation

$$ \partial_t \rho_t = -\bar{\nabla} \cdot j, $$

which is satisfied in the weak sense. This provides a useful heuristic, but as for classical Wasserstein gradient flows [2] the precise, rigorous definition of the tangent space is in Lagrangian form; we note, however, that here we use fluxes instead of velocities. This is not just a superficial difference. Namely, as can be seen in Proposition 2.26, the relation between velocities and fluxes is not linear and thus the velocities do not provide a linear parametrization of the tangent space. We use the argument from [22, Theorem 5.21] to characterize the tangent space $T_\rho P_2(\mathbb{R}^d)$ in more detail.

**Proposition 2.26.** (Tangent fluxes have almost gradient velocities) Let $\mu \in \mathcal{M}^+(\mathbb{R}^d)$ satisfy Assumptions (A1) and (A2), and $\rho \in P_2(\mathbb{R}^d)$. Then, it holds that $j \in T_\rho P_2(\mathbb{R}^d)$ if and only if $j \in \mathcal{M}(G)$ with $j^+ \ll \gamma_1$, $j^- \ll \gamma_2$, and $v^+ := \frac{dj^+}{d\gamma_1}$, $v^- := \frac{dj^-}{d\gamma_2}$ satisfy, for $v := v^+ - v^- : G \to \mathbb{R}$, the relation

$$ v \in \overline{\{\bar{\nabla}\varphi : \varphi \in C_c^\infty(\mathbb{R}^d)\}}^{L^2(\eta, \hat{\gamma}^v)}, \quad \text{where} \quad d\hat{\gamma}^v = \chi_{\{v>0\}} d\gamma_1 + \chi_{\{v<0\}} d\gamma_2. \qquad (2.28) $$

**Proof.** If $\mathcal{A}(\mu; \rho, j) < \infty$, then by Lemma 2.6 it holds for some $v \in \mathcal{V}^{\text{as}}(G)$ that

$$ \begin{aligned} d\tilde{j}(x, y) &= v(x, y)_+ d\gamma_1(x, y) - v(x, y)_- d\gamma_1(y, x) \\ &= v(x, y) d\gamma_+(x, y) - v(y, x) d\gamma_+(y, x), \end{aligned} $$

where $\gamma_+ = \gamma_1|_{J^+}$, with $J^+ = \text{supp } j^+$, and we used that $(J^+)^\top = \text{supp } j^-$. Then, by recalling the definition of the norm on $L^2(\eta \gamma_1)$ from (2.1),

$$ \mathcal{A}(\mu; \rho, j) = 2 \|v_+\|_{L^2(\eta \gamma_1)}^2 = 2 \|v\|_{L^2(\eta \gamma_+)}^2. $$
By using the relation between $j$ and $v$ from above, we can rewrite the divergence $\bar{\nabla} \cdot j$ in weak form for any $\psi \in C_c^\infty(\mathbb{R}^d)$:

$$ \frac{1}{2} \iint_G \bar{\nabla} \psi \, \eta \, d\mathbf{j} = \iint_G \bar{\nabla} \psi \, v_+ \, \eta \, d\gamma_1 = \iint_G \bar{\nabla} \psi \, v \, \eta \, d\gamma_+. $$

Now, the characterization (2.27) of $j \in T_\rho \mathcal{P}_2(\mathbb{R}^d)$ is equivalent to

$$ \iint_G |v|^2 \eta \, d\gamma_+ \le \iint_G |v+w|^2 \eta \, d\gamma_+ \quad \text{for all } w \in \mathcal{V}^{\text{as}}(G) \text{ so that} \\ \iint_G \bar{\nabla} \psi \, w \, \eta \, d\gamma_+ = 0 \quad \forall \psi \in C_c^\infty(\mathbb{R}^d). $$

Hence $v^+$ belongs to the closure of $\{\overline{\nabla}\varphi : \varphi \in C_c^\infty(\mathbb{R}^d)\}$ in $L^2(\eta; \gamma_+)$. From the antisymmetry of $v$ follows that $v^-$ belongs to the closure of $\{\overline{\nabla}\varphi : \varphi \in C_c^\infty(\mathbb{R}^d)\}$ in $L^2(\eta; \gamma_-)$. Thus, the conclusion follows from the identity $\gamma_+ + \gamma_+^\top = \hat{\gamma}^v$ on $G$. $\square$

**Remark 2.27.** Proposition 2.26 shows that for $\mu$ as in its statement, $\rho \in \mathcal{P}_2(\mathbb{R}^d)$ and $j$ chosen from a dense subset of $T_\rho \mathcal{P}_2(\mathbb{R}^d)$, there exists a measurable $\varphi : \mathbb{R}^d \to \mathbb{R}$ such that we have the identity

$$ \mathcal{A}(\mu; \rho, j) = \mathcal{A}(\mu; \rho, \bar{\nabla}\varphi \gamma_1) = \iint_G |(\bar{\nabla}\varphi)_+|^2 \eta \, d\gamma_1. $$

Finally, we provide an interesting property of absolutely continuous curves.

**Proposition 2.28.** (Absolutely continuous curves stay supported on $\mu$) Let $\mu \in \mathcal{M}^+(\mathbb{R}^d)$ satisfy Assumptions (A1) and (A2) and $\rho \in \text{AC}([0, T], (\mathcal{P}_2(\mathbb{R}^d), T_\rho))$ be such that $\text{supp} \rho_0 \subseteq \text{supp} \mu$. Then, for all $t \in [0, T]$, it holds $\text{supp} \rho_t \subseteq \text{supp} \mu$.

**Proof.** Since $(\rho_t)_{t \in [0, T]}$ is absolutely continuous, there exists by Proposition 2.25 a unique family $(j_t)_{t \in [0, T]}$ such that $(\rho, j) \in \text{CE}_T$ and $j_t \in T_{\rho_t} \mathcal{P}_2(\mathbb{R}^d) \subseteq \mathcal{M}_{\gamma_{1,t}}^{\text{as}}(G)$, where $\gamma_{1,t} = \rho_t \otimes \mu$, and $|\rho_t|^2 = \mathcal{A}(\mu; \rho_t, j_t)$ for a.e. $t \in [0, T]$. In particular, by Lemma 2.6, there exists a measurable family $(v_t)_{t \in [0, T]} \subset \mathcal{V}^{\text{as}}(G)$ such that

$$ d\mathbf{j}_t(x, y) = v_t(x, y)_+ d\rho_t(x) d\mu(y) - v_t(x, y)_- d\mu(x) d\rho_t(y). $$

Without loss of generality, let $(\rho_t)_{t \in [0, T]}$ be the weakly continuous curve from Lemma 2.15 satisfying, for any test function $\varphi \in C_c^\infty(\mathbb{R}^d)$ and $t \in [0, T]$,

$$ \begin{aligned} \int_{\mathbb{R}^d} \varphi(x) \, d\rho_t(x) &= \int_{\mathbb{R}^d} \varphi(x) \, d\rho_0(x) + \frac{1}{2} \int_0^t \iint_G \bar{\nabla}\varphi(x, y)\eta(x, y) \, d\mathbf{j}_s(x, y) \, ds \\ &= \int_{\mathbb{R}^d} \varphi(x) \, d\rho_0(x) \\ &\quad + \int_0^t \iint_G \bar{\nabla}\varphi(x, y)v_s(x, y)_+\eta(x, y) \, d(\rho_s \otimes \mu)(x, y) \, ds. \end{aligned} $$
Now, let $\varphi \in C_c^\infty(\mathbb{R}^d)$ with $\varphi \ge 0$ and $\text{supp} \varphi \subseteq \mathbb{R}^d \setminus \text{supp} \mu$. Then, for all $t \in [0, T]$, it holds

$$ \frac{1}{2} \int_{\mathbb{R}^d} \varphi(x) \, d\rho_t(x) = - \int_0^t \iint_G \varphi(x) v_s(x, y)_+ \eta(x, y) \, d(\rho_s \otimes \mu)(x, y) \, ds \le 0, $$

which implies that $\text{supp} \rho_t \subseteq \text{supp} \mu$, since $\rho_t \in \mathcal{P}(\mathbb{R}^d)$ is in particular a non-negative measure for all $t \in [0, T]$ by Lemma 2.15. $\square$

### 3. Nonlocal Nonlocal-Interaction Equation

In this section we consider gradient flows in the spaces of probability measures $\mathcal{P}_2(\mathbb{R}^d)$ endowed with the nonlocal transportation quasi-metric $\mathcal{T}_\mu$, defined by (2.22). From now until Section 3.4 (excluded) we fix $\mu \in \mathcal{M}^+(\mathbb{R}^d)$ satisfying (A1) and (A2), unless otherwise specified. For this reason we shall use the simplifications $\mathcal{A}(\rho, j)$ for $\mathcal{A}(\mu; \rho, j)$ and $\mathcal{T}$ for $\mathcal{T}_\mu$.

In this section investigate the nonlocal nonlocal-interaction equation (NL$^2$IE) as a gradient flow with respect to the metric $\mathcal{T}$. We restate it in a one-line form and note that from now on we consider the external potential $P \equiv 0$. The extension to $P \ne 0$ is straightforward; see Remark 3.2. Thus,

$$ \partial_t \rho_t(x) + \int_{\mathbb{R}^d} \bar{\nabla}(K * \rho)(x, y) - \eta(x, y) \rho_t(x) \, d\mu(y) \quad (\text{NL}^2\text{IE}) \\ - \int_{\mathbb{R}^d} \bar{\nabla}(K * \rho)(x, y) + \eta(x, y) \, d\rho_t(y) = 0. $$

In the classical setting of gradient flows in the spaces of probability measures endowed with the Wasserstein metric [2, 10], the nonlocal-interaction equation

$$ \partial_t \rho_t + \nabla \cdot (\rho_t \nabla (K * \rho_t)) = 0 \quad (3.1) $$

is the gradient flow of the *nonlocal-interaction energy*

$$ \mathcal{E}(\rho) = \frac{1}{2} \iint_{\mathbb{R}^d \times \mathbb{R}^d} K(x, y) \, d\rho(x) \, d\rho(y). \quad (3.2) $$

We start by discussing the geometry of (NL$^2$IE) and interpret it as the gradient flow of (3.2) in the infinite-dimensional Finsler manifold of measures endowed with the Finsler metric associated to $\mathcal{T}$. Following this, we develop a framework of gradient flows in the quasi-metric space $\mathcal{T}$, which extends the setup of gradient flows in metric spaces [2] to quasi-metric spaces. In particular, we build the existence theory for (NL$^2$IE) based on this approach.

Above, for simplicity, (NL$^2$IE) was written for $\rho \ll \mu$, where we recall that we used the notation $\rho$ to denote both the measure and the density with respect to $\mu$. Our framework, however, also applies to the case when $\rho$ is not absolutely continuous with respect to $\mu$. The general weak form of (NL$^2$IE) is obtained in terms of the nonlocal continuity equation as introduced in Section 2.3. Specifically, we have
**Definition 3.1.** A curve $\rho: [0, T] \to \mathcal{P}_2(\mathbb{R}^d)$ is called a weak solution to (NL$^2$IE) if, for the flux $j: [0, T] \to \mathcal{M}(G)$ defined by

$$d\boldsymbol{j}_t(x, y) = \bar{\nabla} \frac{\delta \boldsymbol{\varepsilon}}{\delta \rho}(x, y) - d\rho_t(x) d\mu(y) - \bar{\nabla} \frac{\delta \boldsymbol{\varepsilon}}{\delta \rho}(x, y)_+ d\rho_t(y) d\mu(x),$$

the pair $(\rho, j)$ is a weak solution to the continuity equation

$$\partial_t \rho_t + \bar{\nabla} \cdot \boldsymbol{j}_t = 0 \quad \text{on } [0, T] \times \mathbb{R}^d,$$

according to Definition 2.14.

Here we list the assumptions on the interaction kernel $K: \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}$ we refer to throughout this section:

(K1) $K \in C(\mathbb{R}^d \times \mathbb{R}^d)$;

(K2) $K$ is symmetric, i.e., $K(x, y) = K(y, x)$ for all $(x, y) \in \mathbb{R}^d \times \mathbb{R}^d$;

(K3) $K$ is $L$-Lipschitz near the diagonal and at most quadratic far away, that is there exists some $L \in (0, \infty)$ such that, for all $(x, y), (x', y') \in \mathbb{R}^d \times \mathbb{R}^d$,

$$|K(x, y) - K(x', y')| \le L \left( |(x, y) - (x', y')| \lor |(x, y) - (x', y')|^2 \right).$$

**Remark 3.2.** Assumption (K3) implies that, for some $C > 0$ and all $x, y \in \mathbb{R}^d$,

$$|K(x, y)| \le C (1 + |x|^2 + |y|^2); \qquad (3.3)$$

indeed, for fixed $(x', y') \in \mathbb{R}^d \times \mathbb{R}^d$, (K3) yields

$$|K(x, y)| - |K(x', y')| \le L (1 \lor 2 (|(x, y)|^2 + |(x', y')|^2)),$$

and bounding the maximum $(\lor)$ by the sum, we arrive at $|K(x, y)| \le L + 2L (|(x', y')|^2 + |(x, y)|^2) + |K(x', y')|$, which gives (3.3) with $C = 2L(1 + |(x', y')|^2) + |K(x', y')|$. We notice, by the way, that the bound (3.3) implies that $\mathcal{E}: \mathcal{P}_2(\mathbb{R}^d) \to \mathbb{R}$ is proper with domain equal to $\mathcal{P}_2(\mathbb{R}^d)$.

As mentioned previously, the theory in this section can be easily extended to energies of the form (1.5) including potential energies $\mathcal{E}_P(\rho) = \int_{\mathbb{R}^d} P d\rho$ for some external potential $P: \mathbb{R}^d \to \mathbb{R}$ satisfying a local Lipschitz condition with at-most-quadratic growth at infinity; that is, similarly to (K3), there exists $L \in (0, \infty)$ so that for all $x, y \in \mathbb{R}^d$ we have

$$|P(x) - P(y)| \le L (|x - y| \lor |x - y|^2).$$

We now show that, under the above assumptions on the interaction potential $K$, we have narrow continuity of the energy.

**Proposition 3.3.** (Continuity of the energy) Let the interaction potential $K$ satisfy Assumptions (K1)–(K3). Then, for any sequence $(\rho^n)_n \subset \mathcal{P}_2(\mathbb{R}^d)$ such that $\rho^n \to \rho$ as $n \to \infty$ for some $\rho \in \mathcal{P}_2(\mathbb{R}^d)$, we have

$$\lim_{n \to \infty} \mathcal{E}(\rho^n) = \mathcal{E}(\rho).$$
**Proof.** Let $(\rho^n)_n \subset \mathcal{P}_2(\mathbb{R}^d)$ and $\rho \in \mathcal{P}_2(\mathbb{R}^d)$ be such that $\rho^n \to \rho$ as $n \to \infty$. For all $R > 0$, we write $\bar{B}_R$ the closed ball of radius $R$ centered at the origin in $(\mathbb{R}^d)^2$ and $\varphi_R: (\mathbb{R}^d)^2 \to \mathbb{R}$ a continuous function such that $\varphi_R(z) = 1$ for all $z \in \bar{B}_R$, $\varphi_R(z) = 0$ for all $z \in (\mathbb{R}^d)^2 \setminus \bar{B}_{2R}$, and $\varphi_R(z) \le 1$ for all $z \in (\mathbb{R}^d)^2$. For all $R > 0$, we then set $K_R = \varphi_R K$ and

$$ \mathcal{E}_R(v) = \frac{1}{2} \iint_{\mathbb{R}^d \times \mathbb{R}^d} K_R(x, y) \, dv(y) \, dv(x) \quad \text{for all } v \in \mathcal{P}_2(\mathbb{R}^d). $$

Since $(\rho^n)_n$ converges narrowly to $\rho$ as $n \to \infty$ and $K_R$ is bounded and continuous, we get

$$ \mathcal{E}_R(\rho^n) \to \mathcal{E}_R(\rho) \quad \text{as } n \to \infty. $$

Furthermore, since $K_R \to K$ pointwise as $R \to \infty$, $|K_R| \le |K|$ for all $R > 0$, the domain of $\mathcal{E}$ is $\mathcal{P}_2(\mathbb{R}^d)$ and $\rho \in \mathcal{P}_2(\mathbb{R}^d)$, we also have

$$ \mathcal{E}_R(\rho) \to \mathcal{E}(\rho) \quad \text{as } R \to \infty $$

by the Lebesgue dominated convergence theorem. Similarly, we also have

$$ \mathcal{E}_R(\rho^n) \to \mathcal{E}(\rho^n) \quad \text{as } R \to \infty \text{ for all } n \in \mathbb{N}. $$

By a diagonal argument, we deduce the result. $\square$

### 3.1. Identification of the Gradient in Finsler Geometry

Since the nonlocal upwind transportation cost $\mathcal{T}$ is only a quasi-metric, the underlying structure of $\mathcal{P}_2(\mathbb{R}^d)$ does not have the formal Riemannian structure as it does in the classical gradient flow theory, but a Finslerian structure instead. This highlights the fact that at every point $\rho \in \mathcal{P}_2(\mathbb{R}^d)$ the tangent space $T_\rho \mathcal{P}_2(\mathbb{R}^d)$ is not a Euclidean space, but rather a manifold in its own right.

In this section we provide calculations, in the spirit of Otto's calculus, that characterize the gradient descent in the infinite-dimensional Finsler manifold of probability measures endowed with the nonlocal transportation quasi-metric $\mathcal{T}$. To keep the following considerations simple, we assume that $\rho$ is a given probability measure which is absolutely continuous with respect to $\mu$. In this way, we avoid the need to introduce yet another measure $\lambda \in \mathcal{M}^+(G)$ with respect to which all of the occurring measures are absolutely continuous, similar to how we proceeded in Definition 2.3 for the action. This restriction is done solely to make the presentation clearer and highlight the geometric structure. Hence any flux $j$ of interest is absolutely continuous with respect to $\mu \otimes \mu$ and we can think of $j$ via its density with respect to $\mu \otimes \mu$, which we shall denote by $j$ (using a letter which is not bold).

At every tangent flux $j \in T_\rho \mathcal{P}_2(\mathbb{R}^d)$ we define an inner product $g_{\rho,j}: T_\rho \mathcal{P}_2(\mathbb{R}^d) \times T_\rho \mathcal{P}_2(\mathbb{R}^d) \to \mathbb{R}$ by

$$ g_{\rho,j}(j_1, j_2) = \frac{1}{2} \iint_G j_1(x, y) j_2(x, y) \eta(x, y) $$
$$ \times \left( \frac{\chi_{[j>0]}(x, y)}{\rho(x)} + \frac{\chi_{[j<0]}(x, y)}{\rho(y)} \right) d\mu(x) d\mu(y), \quad (3.4) $$

where {$j > 0$} is an abbreviation for {$ (x, y) \in G : j(x, y) > 0 $} and similarly for {$j < 0$}. The ratios are well-defined since $\rho$ cannot be zero where $j$ is not zero. We note that this is the bilinear form that corresponds to the quadratic form defining the action (see Definition 2.3 and Remark 2.5); namely,

$$ g_{\rho,j}(j, j) = \mathcal{A}(\mu; \rho, j). $$

We refer the reader to “Appendix A” for a derivation of this inner product from a Minkowski norm on $T_{\rho}P_2(\mathbb{R}^d)$ as it is required in Finsler geometry. We recall that from Proposition 2.26 a dense subset of tangent-fluxes $j$ are characterized by the existence of a potential $\varphi \in C_c^{\infty}(\mathbb{R}^d)$ such that, for $\mu \otimes \mu$-a.e. $(x, y) \in G$,

$$ j(x, y) = \bar{\nabla}\varphi(x, y) \left( \rho(x)\chi_{\{\bar{\nabla}\varphi>0\}}(x, y) + \rho(y)\chi_{\{\bar{\nabla}\varphi<0\}}(x, y) \right). \quad (3.5) $$

In this Finsler setting, we now want to determine the direction of steepest descent from $\rho$, for the underlying energy defined in (3.2). The gradient vector of some energy $\mathcal{E}: P(\mathbb{R}^d) \to \mathbb{R}$ at $\rho$, which we denote by grad $\mathcal{E}(\rho)$, is defined as the tangent vector which satisfies

$$ \mathrm{Diff}_{\rho} \mathcal{E}[j] = g_{\rho, \mathrm{grad} \, \mathcal{E}(\rho)}(\mathrm{grad} \, \mathcal{E}(\rho), j) \quad \text{for all } j \in T_{\rho} P_2(\mathbb{R}^d), $$

provided this vector exists and is unique. Here, we use the continuity equation Definition 2.14 to define variations via

$$ \mathrm{Diff}_{\rho} \mathcal{E}[j] = \left. \frac{d}{dt} \right|_{t=0} \mathcal{E}(\tilde{\rho}_t), $$

where $\tilde{\rho}$ is any curve such that $\tilde{\rho}_0 = \rho$ and $\left. \frac{d}{dt} \right|_{t=0} \tilde{\rho}_t = -\bar{\nabla} \cdot j$. From Definition 2.7, due to $\mu \otimes \mu$-absolute continuity of $j$ we have that

$$ -\bar{\nabla} \cdot j(x) = - \int \eta(x, y) j(x, y) d\mu(y) \quad \text{for } \mu\text{-a.e. } x \in \mathbb{R}^d. $$

In the case, when $M$ is a finite-dimensional Finsler manifold, such gradient vector exists and is unique since the mapping $\ell: T_{\rho}M \to (T_{\rho}M)^*$, $j \mapsto g_{\rho,j}(j, \cdot)$, is a bijection; see [18, Proposition 1.9]. For further details into Finsler geometry, we refer the reader to [4,49]. In our case, we can at least claim that the functional $\ell_{\rho}: T_{\rho}P_2(\mathbb{R}^d) \to (T_{\rho}P_2(\mathbb{R}^d))^*$, given for $j \in T_{\rho}P_2(\mathbb{R}^d)$ by

$$ \begin{aligned} j_2 &\mapsto \ell_{\rho}(j)(j_2) = g_{\rho,j}(j, j_2) \\ &= \frac{1}{2} \iint_G j_2(x, y) \eta(x, y) \left( \frac{j(x, y)_+}{\rho(x)} - \frac{j(x, y)_-}{\rho(y)} \right) d\mu(x) d\mu(y), \end{aligned} \quad (3.6) $$

is injective $\eta \mu \otimes \mu$-a.e.; that is, the existence of a gradient implies its uniqueness ($\eta \mu \otimes \mu$-a.e.), in which case we have

$$ \ell_{\rho}(\mathrm{grad} \, \mathcal{E}(\rho)) = \mathrm{Diff}_{\rho} \mathcal{E}. $$
To see the injectivity of (3.6), we first note that $\ell_\rho$ is positively 1-homogeneous by definition. Moreover, we have the following one-sided version of a Cauchy–Schwarz-type estimate

$$ \begin{aligned} \ell_\rho(\mathbf{j})(\mathbf{j}_2) &\le \frac{1}{2} \iint_G \frac{j_2(x, y)_+ j(x, y)_+}{\rho(x)} \eta(x, y) \, d\mu(x) \, d\mu(y) \\ &\quad + \frac{1}{2} \iint_G \frac{j_2(x, y)_- j(x, y)_-}{\rho(y)} \eta(x, y) \, d\mu(x) \, d\mu(y) \\ &\le \sqrt{\ell_\rho(\mathbf{j})(\mathbf{j}) \ell_\rho(\mathbf{j}_2)(\mathbf{j}_2)}. \end{aligned} \quad (3.7) $$

Here, we also used that $\sqrt{ab} + \sqrt{cd} \le \sqrt{(a+c)(b+d)}$ for all $a, b, c, d > 0$. Note that the above inequalities become strict if any of the integrands $j_2(x, y)_+ j(x, y)_-$ or $j_2(x, y)_- j(x, y)_+$ have a contribution. In particular, we could have $\ell_\rho(\mathbf{j})(\mathbf{j}_2) = -\infty$ although the right-hand side is finite. Despite this, we still have equality in (3.7) if and only if $\mathbf{j}_2 = \beta \mathbf{j}_1$ $\eta \mu \otimes \mu$-a.e. for some $\beta \ge 0$.

To prove the injectivity of $\ell_\rho$, let us suppose that $\mathbf{j}_1, \mathbf{j}_2 \in T_\rho \mathcal{P}_2(\mathbb{R}^d)$ are so that $\ell_\rho(\mathbf{j}_1) = \ell_\rho(\mathbf{j}_2)$. If $\mathbf{j}_1 = 0$ or $\mathbf{j}_2 = 0$ $\eta \mu \otimes \mu$-a.e., then $\ell_\rho(\mathbf{j}_1) = \ell_\rho(\mathbf{j}_2)$ implies that $\mathbf{j}_1 = \mathbf{j}_2 = 0$. If both $\mathbf{j}_1$ and $\mathbf{j}_2$ are nonzero, then by the above Cauchy–Schwarz inequality we get

$$ \begin{aligned} 0 < & g_{\rho, j_2}(\mathbf{j}_2, \mathbf{j}_2) = \ell_{\rho}(\mathbf{j}_2)(\mathbf{j}_2) = \ell_{\rho}(\mathbf{j}_1)(\mathbf{j}_2) = g_{\rho, j_1}(\mathbf{j}_1, \mathbf{j}_2) \\ & \le \sqrt{g_{\rho, j_1}(\mathbf{j}_1, \mathbf{j}_1) g_{\rho, j_2}(\mathbf{j}_2, \mathbf{j}_2)}, \end{aligned} $$

which, after dividing by $\sqrt{g_{\rho, j_2}(\mathbf{j}_2, \mathbf{j}_2)}$ yields $g_{\rho, j_2}(\mathbf{j}_2, \mathbf{j}_2) \le g_{\rho, j_1}(\mathbf{j}_1, \mathbf{j}_1)$. Similarly, one gets $g_{\rho, j_1}(\mathbf{j}_1, \mathbf{j}_1) \le g_{\rho, j_2}(\mathbf{j}_2, \mathbf{j}_2)$, from which we get

$$ g_{\rho, j_1}(\mathbf{j}_1, \mathbf{j}_1) = g_{\rho, j_2}(\mathbf{j}_2, \mathbf{j}_2). $$

Hence

$$ \begin{aligned} g_{\rho, j_1}(\mathbf{j}_1, \mathbf{j}_2) &= \ell_{\rho}(\mathbf{j}_1)(\mathbf{j}_2) = \ell_{\rho}(\mathbf{j}_2)(\mathbf{j}_2) = g_{\rho, j_2}(\mathbf{j}_2, \mathbf{j}_2) \\ &= \sqrt{g_{\rho, j_1}(\mathbf{j}_1, \mathbf{j}_1) g_{\rho, j_2}(\mathbf{j}_2, \mathbf{j}_2)}, \end{aligned} $$

which is the equality case in the Cauchy–Schwarz inequality. Therefore, there exists $\beta \ge 0$ such that $\mathbf{j}_2 = \beta \mathbf{j}_1$. By positive 1-homogeneity of $\ell_\rho$ we get $\ell_\rho(\mathbf{j}_2) = \ell_\rho(\beta \mathbf{j}_1) = \beta \ell_\rho(\mathbf{j}_1) = \beta \ell_\rho(\mathbf{j}_2)$, so that $\beta = 1$, since $\ell_\rho(\mathbf{j}_2)(\mathbf{j}_2) \ne 0$. This ends the proof of the claim of injectivity of $\ell_\rho$.

The direction of the steepest descent on Finsler manifolds is in general not $-\mathrm{grad} \, \mathcal{E}(\rho)$, but is defined to be the tangent flux, which we denote by $\mathrm{grad}^{-} \mathcal{E}(\rho)$, such that

$$ -\mathrm{Diff}_{\rho} \mathcal{E}[\mathbf{j}] = g_{\rho, \mathrm{grad}^{-} \mathcal{E}(\rho)} (\mathrm{grad}^{-} \mathcal{E}(\rho), \mathbf{j}) \quad \text{for all } \mathbf{j} \in T_{\rho} \mathcal{P}_{2}(\mathbb{R}^{d}). $$

In other words, we define $\mathrm{grad}^{-} \mathcal{E}(\rho)$ as the tangent vector (provided it exists) such that

$$ \ell_{\rho}(\mathrm{grad}^{-} \mathcal{E}(\rho)) = -\mathrm{Diff}_{\rho} \mathcal{E}. \quad (3.8) $$
Here we clearly see that in general $\mathrm{grad}^{-} \mathcal{E}(\rho) \neq -\mathrm{grad} \mathcal{E}(\rho)$ since $\ell_{\rho}$ is not negatively 1-homogeneous. We can justify that $\mathrm{grad}^{-} \mathcal{E}(\rho)$ indeed corresponds to the direction of steepest descent at $\rho$ via the following criterion, which is analogous to the Riemann case. We first note that if $\mathrm{Diff}_{\rho} \mathcal{E} = 0$ then $\mathrm{grad}^{-} \mathcal{E}(\rho) = 0$. If $\mathrm{Diff}_{\rho} \mathcal{E} \neq 0$ we note that minimizers $j^*$ of

$$j \mapsto \mathrm{Diff}_{\rho} \mathcal{E}[j], \quad \text{with the constraint that } g_{\rho,j}(j, j) = 1,$$

are of the form $j^* = \beta \mathrm{grad}^{-} \mathcal{E}(\rho)$ for some $\beta > 0$. Indeed, using the fact that $\left. \frac{d}{ds} \right|_{s=0} g_{\rho,j+s j_1}(j+s j_1, j+s j_1) = 2g_{\rho,j}(j, j_1)$ for all $j, j_1 \in T_{\rho} \mathcal{P}_2(\mathbb{R}^d)$ [as shown in (A.1) of “Appendix A”] and using the Lagrange multiplier $\beta$ and the functional

$$H(\beta, j) := \mathrm{Diff}_{\rho} \mathcal{E}[j] + \frac{\beta}{2} (g_{\rho,j}(j, j) - 1), \quad j \in T_{\rho} \mathcal{P}_2(\mathbb{R}^d), \quad \beta \in \mathbb{R},$$

yields, for a constrained minimizer $j^*$, the condition

$$\mathrm{Diff}_{\rho} \mathcal{E} = -\beta^* g_{\rho,j^*}(j^*, \cdot) = -\beta^* \ell_{\rho}(j^*). \quad (3.9)$$

By the definition of $j^*$ we have $0 > \mathrm{Diff}_{\rho} \mathcal{E}[j^*] = -\beta^* g_{\rho,j^*}(j^*, j^*)$, which implies that $\beta^* > 0$. By injectivity and positive 1-homogeneity of $\ell_{\rho}$, we get

$$j^* = \ell_{\rho}^{-1} \left( -\frac{1}{\beta^*} \mathrm{Diff}_{\rho} \mathcal{E} \right) = \frac{1}{\beta^*} \ell_{\rho}^{-1} (-\mathrm{Diff}_{\rho} \mathcal{E}) = \frac{1}{\beta^*} \mathrm{grad}^{-} \mathcal{E}(\rho).$$

The gradient flows with respect to $\mathcal{E}$ in the Finsler space $(\mathcal{P}_2(\mathbb{R}^d), \mathcal{T})$ can thus be written

$$\partial_t \rho_t = \bar{\nabla} \cdot \mathrm{grad}^{-} \mathcal{E}(\rho). \quad (3.10)$$

These considerations stay valid for general energy functionals $\mathcal{E}: \mathcal{P}_2(\mathbb{R}^d) \to \mathbb{R}$.

Let us compute the gradient flux for the specific case of the interaction energy (3.2). A direct computation using the symmetry of $K$ and Definition 2.7 gives, for all $j \in T_{\rho} \mathcal{P}_2(\mathbb{R}^d)$,

$$\begin{align*}
- \mathrm{Diff}_{\rho} \mathcal{E}[j] \\
&= \frac{1}{2} \iint_G (-\bar{\nabla}(K * \rho))(x, y) \eta(x, y) j(x, y) \, d\mu(x) \, d\mu(y) \\
&= \frac{1}{2} \iint_G j(x, y) \eta(x, y) \\
&\qquad \times \left( \frac{\rho(x)(-\bar{\nabla}(K * \rho))_+(x, y)}{\rho(x)} - \frac{\rho(y)(-\bar{\nabla}(K * \rho))_-(x, y)}{\rho(y)} \right) \, d\mu(x) \, d\mu(y) \\
&= \frac{1}{2} \iint_G j(x, y) \eta(x, y) (-\bar{\nabla}(K * \rho)(x, y)) \\
&\qquad \times \left( \frac{\rho(x)\chi_{[-\bar{\nabla}K * \rho > 0]}(x, y)}{\rho(x)} + \frac{\rho(y)\chi_{[-\bar{\nabla}K * \rho < 0]}(x, y)}{\rho(y)} \right) \, d\mu(x) \, d\mu(y) \\
&= \ell_{\rho}(\mathrm{grad}^{-} \mathcal{E}(\rho))(j),
\end{align*}$$
where by comparison with (3.6), we observe that grad⁻ *E*(*ρ*) is given for *μ* ⊗ *μ*-a.e. (*x*, *y*) ∈ *G* by

$$
\begin{equation}
\begin{aligned}
\operatorname{grad}^{-} \mathcal{E}(\rho)(x, y) &=-\bar{\nabla}(K * \rho)(x, y) \left(\rho(x) \chi_{[-\bar{\nabla} K * \rho>0]}(x, y)+\rho(y) \chi_{[-\bar{\nabla} K * \rho<0]}(x, y)\right).
\end{aligned}
\tag{3.11}
\end{equation}
$$

This shows by (3.8) the existence and by our previous argument also uniqueness
of grad⁻ *E*(*ρ*). It is easily observed that it has exactly the form (3.5) with the
corresponding potential given by *φ* = −*K* *ρ*.

We conclude this section by mentioning that the Finsler gradient flow structure
of differential equations has been discovered and investigated in other systems; see
[1,41,42].

3.2. Variational Characterization for the Nonlocal Nonlocal-Interaction Equation

Section 3.1 shows that the nonlocal nonlocal-interaction equation (NL²IE) can
in fact be written as the gradient descent of the energy E according to the Finsler
gradient operator; see (3.10) and (3.11). This is why we refer to weak solutions of
(NL²IE) as gradient flows.

In this section we consider ($\mathcal{P}_2(\mathbb{R}^d)$, $\mathcal{T}$) as a quasi-metric space rather than
a Finsler manifold, which allows us to prove rigorous statements more easily. In
particular, we show that the weak solutions of (NL²IE) are curves of maximal
slope for the energy (3.2) in the quasi-metric space ($\mathcal{P}_2(\mathbb{R}^d)$, $\mathcal{T}$) and vice versa.
We then establish the existence and stability of gradient flows using the variational
framework of curves of maximal slope. To develop the variational formulation,
we adapt the approach of [2] to curves of maximal slope in metric spaces to the
quasi-metric space ($\mathcal{P}_2(\mathbb{R}^d)$, $\mathcal{T}$). This requires introducing a *one-sided* version of
the usual concepts from [2] to cope with the asymmetry of the quasi-metric $\mathcal{T}$.

**Definition 3.4.** (One-sided strong upper gradient) A function $h: \mathcal{P}_2(\mathbb{R}^d) \to [0, \infty]$ is a one-sided strong upper gradient for $\mathcal{E}$ if for every $\rho \in \text{AC}([0, T]; (\mathcal{P}_2(\mathbb{R}^d), \mathcal{T}))$ the function $h \circ \rho$ is Borel and

$$
\mathcal{E}(\rho_t) - \mathcal{E}(\rho_s) \geq - \int_s^t h(\rho_\tau) |\rho'_\tau| \, d\tau \quad \text{for all } 0 \leq s \leq t \leq T, \quad (3.12)
$$

where $|\rho'|$ is the metric derivative of $\rho$ as defined in (2.25).

The above one-sided definition is sufficient to characterize the curves of maxi-
mal slope.

**Definition 3.5.** (*Curve of maximal slope*) A curve $\rho \in \text{AC}([0, T]; \mathcal{P}_2(\mathbb{R}^d))$ is a curve of maximal slope for $\mathcal{E}$ with respect to its one-sided strong upper gradient $h$ if and only if $t \mapsto \mathcal{E}(\rho_t)$ is non-increasing and

$$
\mathcal{E}(\rho_t) - \mathcal{E}(\rho_s) + \frac{1}{2} \int_s^t (h(\rho_\tau)^2 + |\rho'_\tau|^2) \, d\tau \le 0 \quad \text{for all } 0 \le s \le t \le T. \tag{3.13}
$$
**Remark 3.6.** Note that by using Young's inequality in (3.12), we get

$$ \mathcal{E}(\rho_t) - \mathcal{E}(\rho_s) + \frac{1}{2} \int_s^t (h(\rho_\tau)^2 + |\rho'_\tau|^2) \, d\tau \geq 0 \quad \text{for all } 0 \leq s \leq t \leq T. $$

Hence, if the curve $(\rho_t)_{t \in [0, T]}$ is a curve of maximal slope for $\mathcal{E}$ with respect to its strong upper gradient $h$, we actually have an equality in (3.13).

Therefore, in order to give a variational characterization of (NL$^2$IE) we need to detect the right one-sided strong upper gradient. As showed in [24], the variation of the energy along the solution to the equation provides the suitable candidate. In what follows we clarify this point as well as the strategy.

We recall that Proposition 2.25 ensures that for any $\rho \in AC([0, T]; (\mathcal{P}_2(\mathbb{R}^d), \mathcal{T}))$ there exists a unique flux $(j_t)_{t \in [0, T]}$ in $T_\rho \mathcal{P}_2(\mathbb{R}^d)$ such that $\int_0^T \mathcal{A}(\rho_t, j_t) \, dt < \infty$, $(\rho, j) \in \text{CE}_T$ and $|\rho'_t|^2 = \mathcal{A}(\rho, j_t)$ for a.e. $t \in [0, T]$. Moreover, according to Lemma 2.6 there exists an antisymmetric measurable vector field $w: [0, T] \times \mathbb{R} \to \mathbb{R}$ such that

$$ d\mathbf{j}_t(x, y) = w_t(x, y)_+ d\gamma_{1,t}(x, y) - w_t(x, y)_- d\gamma_{2,t}(x, y). \quad (3.14) $$

It will be convenient to work directly with this vector field $(w_t)_{t \in [0, T]}$: from now on we write $(\rho, w) \in \text{CE}_T$ for $(\rho, j) \in \text{CE}_T$ as well as $\tilde{\mathcal{A}}(\rho_t, w_t)$ for $\mathcal{A}(\rho_t, j_t)$ according to (2.8). With this convention, we can define a Finsler-type product on velocities in analogy to (3.4) as

$$ \hat{g}_{\rho,w}(u, v) = \frac{1}{2} \iint_G u(x, y) v(x, y) \eta(x, y) \times (\chi_{\{w>0\}}(x, y) d\gamma_1(x, y) + \chi_{\{w<0\}}(x, y) d\gamma_2(x, y)). $$

Note that, under the absolute-continuity assumptions of Section 3.1, by comparing with (3.4) we have that $\hat{g}_{\rho,w}(u, v) = g_{\rho,j}(j_1, j_2)$, where $j_1, j_2$ are obtained from $u, v$ by (3.14), respectively. Moreover, taking (3.6) into account, we also define

$$ \hat{\ell}_{\rho}(w)(v) = \hat{g}_{\rho,w}(w, v). \quad (3.15) $$

Arguing as in (3.7), we arrive at the following one-sided Cauchy–Schwarz inequality:

**Lemma 3.7.** (One-sided Cauchy–Schwarz inequality) For all $v, w \in T_\rho \mathcal{P}_2(\mathbb{R}^d)$ it holds that

$$ \hat{g}_{\rho,w}(w, v) \leq \sqrt{\hat{g}_{\rho,v}(v, v) \hat{g}_{\rho,w}(w, w)}, \quad (3.16) $$

with equality if and only if, for some $\lambda > 0$, $v(x, y)_+ = \lambda w(x, y)_+$ for $\eta \rho \otimes \mu$-a.e. $(x, y) \in G$ (and thus, by antisymmetry, also $v(x, y)_- = \lambda w(x, y)_-$ for $\eta \mu \otimes \rho$-a.e. $(x, y) \in G$).
**Proof.** Using $v = v_+ - v_-$ and the usual Cauchy–Schwarz inequality in $L^2(\eta \rho \otimes \mu)$, we get

$$ 
\begin{aligned}
\hat{g}_{\rho,w}(w, v) &= \frac{1}{2} \iint_G v(x, y)\eta(x, y) \\
&\quad \times (w(x, y)_+ d\rho(x) d\mu(y) - w(x, y)_- d\mu(x) d\rho(y)) \\
&\le \frac{1}{2} \iint_G v(x, y)_+ w(x, y)_+ \eta(x, y) d\rho(x) d\mu(y) \\
&\quad + \frac{1}{2} \iint_G v(x, y)_- w(x, y)_- \eta(x, y) d\mu(x) d\rho(y) \\
&\le \sqrt{\hat{g}_{\rho,v}(v, v) \hat{g}_{\rho,w}(w, w)}.
\end{aligned}
$$

From the usual Cauchy–Schwarz inequality we have equalities above if and only if there exists $\lambda > 0$ such that $v(x, y)_+ = \lambda w(x, y)_+$ for $\eta\rho \otimes \mu$-a.e. $(x, y) \in G$ and $v(x, y)_- = \lambda w(x, y)_-$ for $\eta\mu \otimes \rho$-a.e. $(x, y) \in G$, since all the contributions are positive. $\square$

Now note that, from the weak formulation of the nonlocal continuity equation (2.15), we have for any $\varphi \in C_c^\infty(\mathbb{R}^d)$ and any $0 \le s < t \le T$ the following chain rule:

$$ 
\begin{aligned}
\int_{\mathbb{R}^d} \varphi(x) \, d\rho_t(x) - \int_{\mathbb{R}^d} \varphi(x) \, d\rho_s(x) \\
&= \frac{1}{2} \int_s^t \iint_G \bar{\nabla} \varphi(x, y) \, \eta(x, y) \, d\mathbf{j}_\tau(x, y) \, d\tau \\
&= \frac{1}{2} \int_s^t \iint_G \bar{\nabla} \varphi(x, y) \, \eta(x, y) \\
&\quad \times (w_\tau(x, y)_+ d\gamma_{1,\tau}(x, y) - w_\tau(x, y)_- d\gamma_{2,\tau}(x, y)) \, d\tau \\
&= \frac{1}{2} \int_s^t \iint_G \bar{\nabla} \varphi(x, y) w_\tau(x, y) \, \eta(x, y) \\
&\quad \times (\chi_{\{w>0\}} d\gamma_{1,\tau}(x, y) + \chi_{\{w<0\}} d\gamma_{2,\tau}(x, y)) \, d\tau \\
&= \int_s^t \hat{g}_{\rho_\tau, w_\tau}(w_\tau, \bar{\nabla} \varphi) \, d\tau = \int_s^t \hat{\ell}_\rho(w_\tau)(\bar{\nabla} \varphi) \, d\tau. \quad (3.17)
\end{aligned}
$$

Moreover, we still have the identification of the product $\hat{g}$ with the action in the form of Lemma 2.6,

$$ 
\begin{aligned}
\hat{g}_{\rho_t, w_t}(w_t, w_t) &= \frac{1}{2} \iint_G w_t(x, y)^2 \eta(x, y) \\
&\quad \times (\chi_{\{w>0\}}(x, y) d\gamma_{1,t}(x, y) + \chi_{\{w<0\}}(x, y) d\gamma_{2,t}(x, y)) \\
&= \frac{1}{2} \iint_G w_t(x, y)_+^2 \eta(x, y) d\gamma_{1,t}(x, y) \\
&\quad + \frac{1}{2} \iint_G w_t(x, y)_-^2 \eta(x, y) d\gamma_{2,t}(x, y)
\end{aligned}
$$
$$ 
\begin{aligned}
&= \frac{1}{2} \iint_G \left( w_t(x, y)_+^2 + w_t(y, x)_-^2 \right) \eta(x, y) d\gamma_{1,t}(x, y) \\
&= \hat{\mathcal{A}}(\rho_t, w_t),
\end{aligned} 
\quad (3.18) $$

which shows that the action is the norm with respect to the Finsler structure.

A crucial step toward the variational characterization of (NL²IE) mentioned above is to obtain the chain rule (3.17) for the energy functional (3.2), which is done in Proposition 3.10 below by a suitable regularization. As a consequence, by using the one-sided Cauchy–Schwarz inequality from Lemma 3.7, we obtain in Corollary 3.11 that the square root $\sqrt{D}$ of the local slope, defined below in (3.19), is a one-sided strong upper gradient for $\mathcal{E}$ with respect to the quasi-metric $T$ in the sense of Definition 3.4, where $|\rho_t'|^2 = \hat{\mathcal{A}}(\rho_t, w_t) = \hat{g}_{\rho_t, w_t}(w_t, w_t)$ for a.e. $t \in [0, T]$ due to Proposition 2.25 and (3.18). This allows us to define the De Giorgi functional, which provides the characterization of weak solutions as curves of maximal slope.

**Definition 3.8.** (Local slope and De Giorgi functional) For any $\rho \in \mathcal{P}_2(\mathbb{R}^d)$, let the local slope at $\rho$ be given by

$$ \mathcal{D}(\rho) := \hat{g}_{\rho, -\overline{\nabla} \frac{\delta \mathcal{E}}{\delta \rho}} \left( -\overline{\nabla} \frac{\delta \mathcal{E}}{\delta \rho}, -\overline{\nabla} \frac{\delta \mathcal{E}}{\delta \rho} \right). \quad (3.19) $$

For any $\rho \in \text{AC}([0, T]; (\mathcal{P}_2(\mathbb{R}^d), \mathcal{T}))$, the De Giorgi functional at $\rho$ is defined as

$$ \mathcal{G}_T(\rho) := \mathcal{E}(\rho_T) - \mathcal{E}(\rho_0) + \frac{1}{2} \int_0^T (\mathcal{D}(\rho_\tau) + |\rho'_\tau|^2) \, d\tau. \quad (3.20) $$

When the dependence on the base measure $\mu$ needs to be explicit, the local slope and the De Giorgi functional are denoted by $\mathcal{D}(\mu; \rho)$ and $\mathcal{G}_T(\mu; \rho)$, respectively.

If the potential $K$ satisfies Assumptions (K1)–(K3), we note that whenever $\rho$ is a weak solution to (NL²IE) and $\rho \in \text{AC}([0, T]; \mathcal{P}_2(\mathbb{R}^d))$ the quantity $\mathcal{G}_T(\rho)$ is finite; indeed, the domain of the energy is all of $\mathcal{P}_2(\mathbb{R}^d)$ and Proposition 2.25 yields that both the local slope (since it is equal to the action of $(\rho, j)$, where $j$ is given in Definition 3.1) and metric derivative are finite.

We are ready to state our main theorem.

**Theorem 3.9.** Suppose that $\mu$ satisfies Assumptions (A1) and (A2) and $K$ satisfies Assumptions (K1)–(K3). A curve $(\rho_t)_{t \in [0, T]} \subset \mathcal{P}_2(\mathbb{R}^d)$ is a weak solution to (NL²IE) according to Definition 3.1 if and only if $\rho$ belongs to $\text{AC}([0, T]; (\mathcal{P}_2(\mathbb{R}^d), \mathcal{T}))$ and is a curve of maximal slope for $\mathcal{E}$ with respect to $\sqrt{D}$ in the sense of Definition 3.5, that is, satisfies

$$ \mathcal{G}_T(\rho) = 0, \quad (3.21) $$

where $\mathcal{G}_T$ is the De Giorgi functional as given in Definition 3.8.
Note that in the above theorem, the implicit assumption that $\sqrt{D}$ is a one-sided strong upper gradient for $\mathcal{E}$ is made; this is in fact true thanks to Corollary 3.11 below. In light of this we can represent the result via the following diagram:

$\rho$ is a weak solution of (NL$^2$IE)

iff $\rho$ is a curve of maximal slope for $\mathcal{E}$ w.r.t. $\sqrt{D}$

iff $G_T(\rho)=0$.

### 3.3. The Chain Rule and Proof of Theorem 3.9

Firstly, we focus on the chain-rule property, which is the main technical step for proving Theorem 3.9.

**Proposition 3.10.** Let $K$ satisfy Assumptions (K1)-(K3). For all $\rho \in AC([0, T]; (P_2(\mathbb{R}^d), T))$ and $0 \le s \le t \le T$ we have the chain-rule identity

$$ \mathcal{E}(\rho_t) - \mathcal{E}(\rho_s) = \int_s^t \hat{g}_{\rho_\tau, w_\tau} \left( w_\tau, \bar{\nabla} \frac{\delta \mathcal{E}}{\delta \rho}(\rho_\tau) \right) d\tau, \quad (3.22) $$

where $(w_t)_{t \in [0, T]}$ is the antisymmetric vector field associated by (2.6) to $(\rho, j) \in CE_T$.

**Proof.** Since the curve $\rho \in AC([0, T]; (P_2(\mathbb{R}^d), T))$, according to Proposition 2.25 there exists a unique family $(j_t)_{t \in [0, T]}$ belonging to $T_\rho P_2(\mathbb{R}^d)$ for a.e. $t \in [0, T]$ such that:

(i) $(\rho, j) \in CE_T$;

(ii) $\int_0^T \sqrt{A(\rho_t, j_t)} dt < \infty$;

(iii) $|\rho'_t|^2 = A(\rho_t, j_t)$ for a.e. $t \in [0, T]$;

(iv) $d\mathbf{j}_t(x, y) = w_t(x, y)_+ d\gamma_{1,t}(x, y) - w_t(x, y)_- d\gamma_{2,t}(x, y)$.

Then the identity (3.22) is equivalent to proving

$$ \mathcal{E}(\rho_t) - \mathcal{E}(\rho_s) = \frac{1}{2} \iint_G \bar{\nabla} \frac{\delta \mathcal{E}}{\delta \rho}(\rho_\tau)(x, y) \eta(x, y) d\mathbf{j}_\tau(x, y) d\tau. \quad (3.23) $$

We proceed by applying two regularization procedures. First, for all $(x, y) \in \mathbb{R}^d \times \mathbb{R}^d$ we define $K^\epsilon(x, y) = K * m_\epsilon(x, y) = \iint_{\mathbb{R}^d \times \mathbb{R}^d} K(z, z')m_\epsilon(x - z, y - z') dz dz'$, where $m_\epsilon(z) = \frac{1}{\epsilon^{2d}}m(\frac{z}{\epsilon})$ for all $z \in \mathbb{R}^{2d}$ and $\epsilon > 0$, where $m$ is a standard mollifier on $\mathbb{R}^{2d}$. We also introduce a smooth cut-off function $\varphi_R$ on $\mathbb{R}^{2d}$ such that $\varphi(z) = 1$ on $B_R$, $\varphi(z) = 0$ on $\mathbb{R}^{2d} \setminus B_{2R}$ and $|\nabla\varphi_R| \le \frac{2}{R}$, where $B_R$ is the ball of radius $R$ in $\mathbb{R}^{2d}$ centered at the origin. We set $K_R^\epsilon := \varphi_R K^\epsilon$ and note that it is a $C_c^\infty(\mathbb{R}^{2d})$ function. We now introduce the approximate energies, indexed by $\epsilon$ and $R$,

$$ \mathcal{E}_R^\epsilon(v) = \frac{1}{2} \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} K_R^\epsilon(x, y) \, dv(y) \, dv(x) \quad \text{for all } v \in P_2(\mathbb{R}^d). $$
Let us extend $\rho$ and $j$ to $[-T, 2T]$ periodically in time, meaning that $\rho_{-s} = \rho_{T-s}$ and $\rho_{T+s} = \rho_s$ for all $s \in (0, T]$ and likewise for $j$. We regularize $\rho$ and $j$ in time by using a standard mollifier $n$ on $\mathbb{R}$ supported on $[-1, 1]$, by setting $n_\sigma(t) = \frac{1}{\sigma^n(t)}$ and

$$ \rho_t^\sigma(A) = n_\sigma * \rho_t(A) = \int_{-\sigma}^\sigma n_\sigma(t-s) \rho_s(A) \, ds, \quad \forall A \subseteq \mathbb{R}^d, $$

$$ j_t^\sigma(U) = n_\sigma * j_t(A) = \int_{-\sigma}^\sigma n_\sigma(t-s) j_s(U) \, ds, \quad \forall U \subset G, $$

for any $\sigma \in (0, T)$; whence $\rho_t^\sigma \in P_2(\mathbb{R}^d)$. Let us now show that the integral of the action is uniformly bounded with respect to $\sigma$. Let $|\lambda| \in M^+(G)$ be such that $\gamma_{1,t}, \gamma_{2,t}, |j_t| \ll |\lambda|$ for all $t \in [0, T]$. Then by using the joint convexity of the function $\alpha$ from (2.5), Jensen's inequality and Fubini's Theorem, we get

$$ \begin{aligned} & \int_0^T A(\rho_t^\sigma, j_t^\sigma) \, dt \\ &= \frac{1}{2} \int_0^T \iint_G \alpha \left( \int_{-\sigma}^\sigma \frac{d\mathbf{j}_{t-s}}{d|\lambda|} n_\sigma(s) \, ds, \int_{-\sigma}^\sigma \frac{d\gamma_{1,t-s}}{d|\lambda|} n_\sigma(s) \, ds \right) \eta \, d|\lambda| \, dt \\ & \quad + \frac{1}{2} \int_0^T \iint_G \alpha \left( -\int_{-\sigma}^\sigma \frac{d\mathbf{j}_{t-s}}{d|\lambda|} n_\sigma(s) \, ds, \int_{-\sigma}^\sigma \frac{d\gamma_{2,t-s}}{d|\lambda|} n_\sigma(s) \, ds \right) \eta \, d|\lambda| \, dt \\ & \le \frac{1}{2} \int_0^T \iint_G \int_{-\sigma}^\sigma \alpha \left( \frac{d\mathbf{j}_{t-s}}{d|\lambda|}, \frac{d\gamma_{1,t-s}}{d|\lambda|} \right) n_\sigma(s) \, ds \, \eta \, d|\lambda| \, dt \\ & \quad + \frac{1}{2} \int_0^T \iint_G \int_{-\sigma}^\sigma \alpha \left( -\frac{d\mathbf{j}_{t-s}}{d|\lambda|}, \frac{d\gamma_{2,t-s}}{d|\lambda|} \right) n_\sigma(s) \, ds \, \eta \, d|\lambda| \, dt \\ &= \int_{-\sigma}^{+\sigma} \int_0^T A(\rho_{t-s}, j_{t-s}) \, dt \, n_\sigma(s) \, ds \\ & \le \int_{-T}^{2T} A(\rho_t, j_t) \, dt = 3 \int_0^T A(\rho_t, j_t) \, dt < \infty. \end{aligned} $$

It is easy to check that $(\rho^\sigma, j^\sigma)$ is still a solution to the nonlocal continuity equation on $[0, T]$. By arguing as in the proof of Proposition 2.17, we get that along subsequences it holds $\rho_t^\sigma \to \tilde{\rho}_t$ as $\sigma \to 0$ for all $t \in [0, T]$ for some curve $(\tilde{\rho}_t)_{t \in [0, T]}$ in $P_2(\mathbb{R}^d)$, and $j^\sigma \to \hat{j}$ in $M_\text{loc}(G \times [0, T])$. With $d\hat{j} := d\tilde{j}_t \, dt$, for some curve $(\tilde{j}_t)_{t \in [0, T]}$ in $M(G)$. Note that $n_\sigma \to \delta_0$ as $\sigma \to 0$, and, as a consequence, $\rho_t^\sigma \to \rho_t$ for all $t \in [0, T]$ in the view of Proposition 2.21. Thus, we actually have $\tilde{\rho} = \rho$ and $\hat{j} = j$ by uniqueness of the limit and the flux, as highlighted above. Using the regularity for $\varepsilon > 0$ and $\sigma > 0$, we get

$$ \begin{aligned} \frac{d}{dt} \mathcal{E}_R^\varepsilon(\rho_t^\sigma) &= \int_{\mathbb{R}^d} (K_R^\varepsilon * \rho_t^\sigma)(x) \partial_t \rho_t^\sigma(x) \, d\mu(x) \\ &= \frac{1}{2} \iint_G \bar{\nabla} (K_R^\varepsilon * \rho_t^\sigma)(x, y) \eta(x, y) \, d\mathbf{j}_t^\sigma(x, y). \end{aligned} $$

For the sake of completeness, we note that the second equality follows from the definition of $CE_T$ by using again a cut-off argument on the function $K_R^\varepsilon * \rho_t^\sigma$. We
omit this step as it is a standard procedure. By integrating in time between $s$ and $t$, with $s \le t$, it follows

$$ 
\begin{aligned}
\mathcal{E}_R^\epsilon(\rho_t^\sigma) - \mathcal{E}_R^\epsilon(\rho_s^\sigma) \\
&= \frac{1}{2} \int_s^t \iint_G \bar{\nabla}(K_R^\epsilon * \rho_\tau^\sigma)(x, y) \eta(x, y) \, d\mathbf{j}_\tau^\sigma(x, y) \, d\tau \\
&= \frac{1}{2} \int_s^t \iint_G \int_{\mathbb{R}^d} (K_R^\epsilon(y, z) - K_R^\epsilon(x, z)) \, d\rho_\tau^\sigma(z) \eta(x, y) \, d\mathbf{j}_\tau^\sigma(x, y) \, d\tau.
\end{aligned} 
\quad (3.24) $$

In order to obtain (3.23) we need to let $\epsilon$ and $\sigma$ go to 0 and $R$ go to $\infty$ in (3.24). The left-hand side is easy to handle since $\rho_t^\sigma \to \rho_t$ as $\sigma \to 0$ for any $t \in [0, T]$, and $K_R^\epsilon \to K_R$ uniformly on compact sets as $\epsilon \to 0$. Finally, by letting $R$ go to $\infty$ we have convergence to $\mathcal{E}(\rho_t)$.

In order to pass to the limit in the right-hand side of (3.24), we use a truncation argument similar to that in the proof of Proposition 2.17. Let $\delta > 0$ and let us set $N_\delta = \bar{B}_{\delta-1} \times \bar{B}_{\delta-1}$, where $B_{\delta-1} = \{x \in \mathbb{R}^d : |x| < \delta^{-1}\}$, and $G_\delta = \{(x, y) \in G : \delta \le |x - y|\}$. We can consider a family $(\varphi_\delta)_{\delta>0} \subset C_c^\infty(\mathbb{R}^d \times G; [0, 1])$ of truncation functions such that, for all $\delta > 0$,

$$ \{\varphi_\delta = 1\} \supseteq \bar{B}_{\delta-1} \times G_\delta \cap N_\delta. $$

Now, we add and subtract $\varphi_\delta$ in the integral on the RHS of (3.24) and we argue as follows. Since $\rho_t^\sigma \otimes \mathbf{j}_t^\sigma \to \rho_t \otimes \mathbf{j}_t$ for any $t \in [0, T]$ as $\sigma \to 0$, and $K_R^\epsilon \to K_R$ uniformly on compact sets as $\epsilon \to 0$, we can pass to the limit in $\sigma$ and $\epsilon$, for any $R$ and $\delta > 0$,

$$ 
\begin{aligned}
& \frac{1}{2} \int_s^t \iint_G \int_{\mathbb{R}^d} \varphi_\delta(z, x, y) (K_R^\epsilon(y, z) - K_R^\epsilon(x, z)) \, d\rho_\tau^\sigma(z) \eta(x, y) \, d\mathbf{j}_\tau^\sigma(x, y) \, d\tau \\
& \to \frac{1}{2} \int_s^t \iint_G \int_{\mathbb{R}^d} \varphi_\delta(z, x, y) (K_R(y, z) - K_R(x, z)) \, d\rho_\tau^\sigma(z) \eta(x, y) \, d\mathbf{j}_\tau^\sigma(x, y) \, d\tau.
\end{aligned} 
\quad (3.25) $$

By using $\varphi_\delta \le 1$, Assumption (K3), Lemma 2.10 with $\Phi(x, y) = |x - y| \vee |x - y|^2$ and (A1), we can bound the modulus of (3.25) for any $\tau \in [s, t]$ by

$$ 
\begin{aligned}
& \frac{1}{2} \iint_G \int_{\mathbb{R}^d} \frac{|K_R(y, z) - K_R(x, z)|}{|x - y| \vee |x - y|^2} \, d\rho_\tau(z) (|x - y| \vee |x - y|^2) \eta(x, y) \, d|\mathbf{j}_\tau|(x, y) \\
& \le L \sqrt{2C_\eta} \mathcal{A}(\rho_t, \mathbf{j}_t).
\end{aligned} 
$$

Hence the integral is uniformly bounded in $\delta$ and $R$, and by the Lebesgue dominated convergence theorem we can pass to the limit in (3.25) in $\delta$ and $R$, obtaining

$$ \frac{1}{2} \int_s^t \iint_G \int_{\mathbb{R}^d} (K(y, z) - K(x, z)) \, d\rho_\tau(z) \eta(x, y) \, d\mathbf{j}_\tau(x, y) \, d\tau. $$
Now, it remains to control the integral involving the term $1 - \varphi_\delta(z, x, y)$ in the integrand. Let us note that, for all $\delta > 0$,

$$ (\mathbb{R}^d \times G) \setminus \{\varphi_\delta = 1\} \subseteq (\bar{B}_{\delta-1}^c \times G) \cup (\mathbb{R}^d \times (G \setminus (G_\delta \cap N_\delta))) =: M_\delta. $$

Using Assumption (K3) and splitting each contribution, we obtain

$$ \begin{aligned} & \left| \iint_G \int_{\mathbb{R}^d} (1 - \varphi_\delta(z, x, y)) (K_R^\epsilon(y, z) - K_R^\epsilon(x, z)) d\rho_t^\sigma(z) \eta(x, y) d\mathbf{j}_t^\sigma(x, y) \right| \\ & \le L \iiint_{M_\delta} (|x-y| \vee |x-y|^2) \eta(x, y) d\mathbf{j}_t^\sigma(x, y) d\rho_t^\sigma(z) \\ & \le L \iiint_{\bar{B}_{\delta-1}^c \times G} (|x-y| \vee |x-y|^2) \eta(x, y) d\mathbf{j}_t^\sigma(x, y) d\rho_t^\sigma(z) \\ & \qquad + 2L \int_{\mathbb{R}^d} d\rho_t^\sigma(z) \iint_{G_\delta^c} (|x-y| \vee |x-y|^2) w_t(x, y) + \eta(x, y) d\rho_t^\sigma(x) d\mu(y) \\ & \qquad + 2L \int_{\mathbb{R}^d} d\rho_t^\sigma(z) \iint_{N_\delta^c} (|x-y| \vee |x-y|^2) w_t(x, y) + \eta(x, y) d\rho_t^\sigma(x) d\mu(y). \end{aligned} $$

Using Lemma 2.10 with $\Phi(x, y) = |x - y| \vee |x - y|^2$, (A1) and the Cauchy–Schwarz inequality with respect to $\eta \rho_t^\sigma \otimes \mu$, the right-hand side in the inequality above can be further bounded by

$$ 4L \sqrt{C_\eta A(\rho_t^\sigma, j_t^\sigma)} \rho_t^\sigma(\bar{B}_{\delta-1}^c) + 2L \sqrt{A(\rho_t^\sigma, j_t^\sigma)} \left( \left( \iint_{G_\delta^c} |x-y|^2 \eta(x,y) d\rho_t^\sigma(x) d\mu(y) \right)^{\frac{1}{2}} + \sqrt{C_\eta \rho_t^\sigma(\bar{B}_{\delta-1}^c)} \right). $$

Thanks to the uniform second moment bound of $\rho_t^\sigma$ from Lemma 2.16 and Assumption (A2), the above terms converge to zero as $\delta \to 0$, which concludes the proof. $\square$

That $\sqrt{D}$ is a one-sided strong upper gradient for $\mathcal{E}$ is an easy consequence of the previous result.

**Corollary 3.11.** For any curve $\rho \in AC([0, T]; (P_2(\mathbb{R}^d), T))$ it holds that

$$ \mathcal{E}(\rho_t) - \mathcal{E}(\rho_s) \geq - \int_s^t \sqrt{D(\rho_\tau)} |\rho'_\tau| \, d\tau \quad \text{for all } 0 \leq s \leq t \leq T, \quad (3.26) $$

i.e., $\sqrt{D}$ is a one-sided strong upper gradient for $\mathcal{E}$ in the sense of Definition 3.4.

**Proof.** Without loss of generality we assume $\int_s^t \sqrt{D(\rho_\tau)} |\rho'_\tau| (\tau) \, d\tau < \infty$, as otherwise the inequality (3.26) is trivially satisfied. We obtain the result as consequence of Proposition 3.10 by applying the one-sided Cauchy–Schwarz inequality (Lemma 3.7) to (3.22) as follows: for any $0 \leq s \leq t \leq T$,

$$ \mathcal{E}(\rho_t) - \mathcal{E}(\rho_s) $$
$$
\begin{align*}
&= \int_s^t \widehat{g}_{\rho_\tau, w_\tau} \left( w_\tau, \bar{\nabla} \frac{\delta \mathcal{E}(\rho_\tau)}{\delta \rho} \right) d\tau = - \int_s^t \widehat{g}_{\rho_\tau, w_\tau} \left( w_\tau, -\bar{\nabla} \frac{\delta \mathcal{E}(\rho_\tau)}{\delta \rho} \right) d\tau \\
&\geq - \int_s^t \sqrt{\widehat{g}_{\rho_\tau, -\bar{\nabla} \frac{\delta \mathcal{E}(\rho_\tau)}{\delta \rho}} \left( -\bar{\nabla} \frac{\delta \mathcal{E}(\rho_\tau)}{\delta \rho}, -\bar{\nabla} \frac{\delta \mathcal{E}(\rho_\tau)}{\delta \rho} \right)} \sqrt{\widehat{g}_{\rho_\tau, w_\tau} (w_\tau, w_\tau)} d\tau \\
&= \int_s^t \sqrt{\mathcal{D}(\rho_\tau)} \sqrt{\hat{\mathcal{A}}(\rho_\tau, w_\tau)} d\tau \\
&= \int_s^t \sqrt{\mathcal{D}(\rho_\tau)} |\rho'|(\tau) d\tau.
\end{align*}
$$

Note that the last two equalities are provided by identity (3.18) and Proposition 2.25.

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

At this point, we have collected all auxiliary results to deduce Theorem 3.9.

**Proof of Theorem 3.9.** Let us start by assuming that $\rho$ is a weak solution to (NL$^2$IE). In view of Definition 3.1, a weak solution is obtained from the weak formulation of the nonlocal continuity equation (2.13) if we set

$$d\mathbf{j}_t(x, y) = \bar{\nabla} \frac{\delta \mathcal{E}}{\delta \rho}(x, y) - d\rho_t(x) d\mu(y) - \bar{\nabla} \frac{\delta \mathcal{E}}{\delta \rho}(x, y) + d\rho_t(y) d\mu(x).$$

Then, by writing $v_t^\mathcal{E}(x, y) = -\bar{\nabla} \frac{\delta \mathcal{E}}{\delta \rho}(x, y)$, it is easy to check that

$$A(\rho_t, \mathbf{j}_t) = \hat{A}(\rho_t, v_t^\mathcal{E}) = \mathcal{D}(\rho_t) < \infty,$$

where the finiteness follows from Assumptions (K3) and (A1), as shown by the computation

$$
\begin{align*}
\mathcal{D}(\rho_t) &= \iint_G |(\bar{\nabla} K * \rho_t(x, y))_-|^2 \eta(x, y) d\rho_t(x) d\mu(y) \\
&\leq \iint_G (\bar{\nabla} K * \rho_t(x, y))^2 \eta(x, y) d\rho_t(x) d\mu(y) \\
&= \iint_G \left( \int_{\mathbb{R}^d} (K(x, z) - K(y, z)) d\rho_t(z) \right)^2 \eta(x, y) d\rho_t(x) d\mu(y) \\
&\leq \iint_G \int_{\mathbb{R}^d} (K(x, z) - K(y, z))^2 d\rho_t(z) \eta(x, y) d\rho_t(x) d\mu(y) \\
&\leq L^2 \int_{\mathbb{R}^d} \iint_G (|x-y|^2 \lor |x-y|^4) \eta(x, y) d\mu(y) d\rho_t(x) d\rho_t(z) \\
&\leq L^2 C_\eta \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} d\rho_t(x) d\rho_t(z) = L^2 C_\eta.
\end{align*}
$$

Thanks to Proposition 2.25, this also proves that $\rho \in AC([0, T], (P_2(\mathbb{R}^d), T))$ and $|\rho'_t|^2 \leq \mathcal{D}(\rho_t)$ for a.e. $t \in [0, T]$. In view of Proposition 3.10, we thus obtain

$$\mathcal{E}(\rho_t) - \mathcal{E}(\rho_s) = \int_s^t \widehat{g}_{\rho_\tau, v_\tau^\mathcal{E}} \left( v_\tau^\mathcal{E}, \bar{\nabla} \frac{\delta \mathcal{E}}{\delta \rho} (\rho_\tau) \right) d\tau$$
$$
\begin{align*}
&= - \int_s^t \widehat{g}_{\rho_\tau, v_\tau^\varepsilon} \left( v_\tau^\varepsilon, -\bar{\nabla} \frac{\delta \mathcal{E}}{\delta \rho} (\rho_\tau) \right) \mathrm{d}\tau \\
&= - \int_s^t \iint_G \left| \bar{\nabla} \frac{\delta \mathcal{E}}{\delta \rho} (x, y) - \right|^2 \eta(x, y) \mathrm{d}\rho_\tau(x) \mathrm{d}\mu(y) \mathrm{d}\tau \\
&= - \int_s^t \mathcal{D}(\rho_\tau) \mathrm{d}\tau \leq - \int_s^t \sqrt{\mathcal{D}(\rho_\tau)} |\rho'_\tau| \mathrm{d}\tau.
\end{align*}
$$

This implies that

(i) the map $t \mapsto \mathcal{E}(\rho_t)$ is non-increasing;

(ii) $\mathcal{E}(\rho_t) - \mathcal{E}(\rho_s) + \frac{1}{2} \int_s^t \mathcal{D}(\rho_\tau) + |\rho'_\tau|^2 \, d\tau = 0$, by Corollary 3.11.

Whence the first part of the theorem follows for $s = 0$ and $t = T$ since $\mathcal{G}_T(\rho) = 0$.

Consider now $\rho \in AC([0, T]; (\mathcal{P}_2(\mathbb{R}^d), \mathcal{T}))$ satisfying the equality (3.21). Let us verify that it is a weak solution of (NL$^2$IE). By Proposition 2.25 there exists a unique family $(j_t)_{t \in [0, T]}$ in $T_{\rho_t} \mathcal{P}_2(\mathbb{R}^d)$ such that $(\rho, j) \in \text{CE}_T, \int_0^T \sqrt{\mathcal{A}(\rho_t, j_t)} \, dt < \infty$ and $|\rho'_t|^2 = \mathcal{A}(\rho_t, j_t)$ for a.e. $t \in [0, T]$. Moreover, by Lemma 2.6 we find an antisymmetric measurable vector field $w: [0, T] \times G \to \mathbb{R}$ such that

$$
dj_t(x, y) = w_t(x, y)_+ d\gamma_{1,t}(x, y) - w_t(x, y)_- d\gamma_{2,t}(x, y).
$$

Thanks to Proposition 3.10, by applying the one-sided Cauchy–Schwarz, using the identification (3.18), the definition of the local slope (3.19) and Young inequality, we get

$$
\begin{align*}
\mathcal{E}(\rho_T) - \mathcal{E}(\rho_0) &= \int_0^T \widehat{g}_{\rho_\tau, w_\tau} \left( w_\tau, \bar{\nabla} \frac{\delta \mathcal{E}}{\delta \rho} (\rho_\tau) \right) \mathrm{d}\tau \\
&= - \int_0^T \widehat{g}_{\rho_\tau, w_\tau} \left( w_\tau, -\bar{\nabla} \frac{\delta \mathcal{E}}{\delta \rho} (\rho_\tau) \right) \mathrm{d}\tau \\
&\geq - \int_0^T \sqrt{\mathcal{D}(\rho_\tau)} \sqrt{\mathcal{A}(\rho_\tau, j_\tau)} \mathrm{d}\tau = - \int_0^T \sqrt{\mathcal{D}(\rho_\tau)} |\rho'_\tau| \mathrm{d}\tau \\
&\geq - \frac{1}{2} \int_0^T \mathcal{D}(\rho_\tau) \mathrm{d}\tau - \frac{1}{2} \int_0^T |\rho'_\tau|^2 \mathrm{d}\tau.
\end{align*}
$$

Thanks to the equality (3.21), we actually have that the above inequalities are equalities, which holds if and only if $w_t(x, y) = -\bar{\nabla} \frac{\delta \mathcal{E}}{\delta \rho}(x, y)$ for a.e. $t \in [0, T]$ and $\gamma_{1,t}$-a.e. $(x, y) \in G$. Hence $(\rho, j) \in \text{CE}_T$ with $w = -\bar{\nabla} \frac{\delta \mathcal{E}}{\delta \rho}$, that is, $\rho$ is a weak solution to (NL$^2$IE). $\square$

3.4. Stability and Existence of Weak Solutions

Theorem 3.9 provides a characterization of (weak) solutions to (NL$^2$IE) as minimizers of $\mathcal{G}_T$ attaining the value 0. The direct method of calculus of variations gives existence of minimizers of $\mathcal{G}_T$. However, it is not clear a priori whether they attain the value 0 and are thus actually weak solutions to (NL$^2$IE). Hence we prove compactness and stability of gradient flows (see Theorem 3.14) and approximate
the desired problem by discrete problems for which the existence of solutions is easy to show; see the proof of Theorem 3.15. We start by proving that the local slope $\mathcal{D}$ is narrowly lower semicontinuous jointly in its arguments, $\mu$ and $\rho$; see Lemma 3.12. We then establish the compactness coming from a uniform control of the De Giorgi functional $\mathcal{G}_T$, as well as its joint narrow lower semicontinuity (see Lemma 3.13), which we prove using compactness in $CE_T$ and the joint narrow lower semicontinuity of the action (see Proposition 2.17) and of the local slope. (See also [48, Theorem 2] for an analogous strategy.)

In Theorem 3.14 we prove one of our main results, namely that the functional $\mathcal{G}_T$ is stable under variations in base measures, defining the vertices of the graph, and absolutely continuous curves. A particular consequence of this theorem is that weak solutions to (NL$^2$IE) with respect to graphs defined by random samples of a measure $\mu$ converge to weak solutions to (NL$^2$IE) with respect to $\mu$; see Remark 3.17.

The existence of weak solutions of (NL$^2$IE) (and thus gradient flows) with respect to $\mathcal{E}$ proved in Theorem 3.15 shows that, indeed, the De Giorgi functional (3.20) corresponding to an interaction potential $K$ satisfying (K1)–(K3) admits a minimizer when $\mu(\mathbb{R}^d)$ is finite.

**Lemma 3.12.** Let $(\mu^n)_n \subset \mathcal{M}^+(\mathbb{R}^d)$ and suppose that $(\mu^n)_n$ narrowly converges to $\mu$. Assume that the base measures $(\mu^n)_n$ and $\mu$ are such that (A1) and (A2) hold uniformly in $n$, and let $K$ satisfy Assumptions (K1)–(K3). Let moreover $(\rho^n)_n$ be a sequence such that $\rho^n \in \mathcal{P}_2(\mathbb{R}^d)$ for all $n \in \mathbb{N}$ and $\rho^n \to \rho$ as $n \to \infty$ for some $\rho \in \mathcal{P}_2(\mathbb{R}^d)$. Then

$$ \liminf_{n \to \infty} \mathcal{D}(\mu^n; \rho^n) \geq \mathcal{D}(\mu; \rho). $$

**Proof.** For every $n \in \mathbb{N}$ we set $u^n = \bar{\nabla} K * \rho^n$. Furthermore, we write $u = \bar{\nabla} K * \rho$ and define $g: \mathbb{R} \to \mathbb{R}$ by $g(x) = (x_+)^2$ for all $x \in \mathbb{R}$. Then note that $g$ is convex and continuous, and

$$ \mathcal{D}(\mu^n; \rho^n) = \iint_G g(u^n(x, y))\eta(x, y) \, d\rho^n(x) \, d\mu^n(y), $$

and, similarly,

$$ \mathcal{D}(\mu; \rho) = \iint_G g(u(x, y))\eta(x, y) \, d\rho(x) \, d\mu(y). $$

We want to use [2, Theorem 5.4.4 (ii)] to prove the desired $\liminf$ inequality. Observe that $u^n \in L^2(\eta \gamma_1^n)$ and $u \in L^2(\eta \gamma_1)$; indeed, (K3) and (A1) give

$$ \begin{align*} \iint_G u^n(x, y)^2 \eta(x, y) \, d\gamma_1^n(x, y) & \\ &= \iint_G (K * \rho^n(y) - K * \rho^n(x))^2 \eta(x, y) \, d\gamma_1^n(x, y) \\ &\leq L^2 C_\eta, \end{align*} $$
and, similarly, for $u$. Let now $\varphi \in C_c^\infty(G)$. We have

$$
\begin{align*}
& \iint_G u^n(x, y) \varphi(x, y) \eta(x, y) \, d\gamma_1^n(x, y) \\
&= \iint_G \left( \int_{\mathbb{R}^d} K(y, z) \, d\rho^n(z) - \int_{\mathbb{R}^d} K(x, z) \, d\rho^n(z) \right) \varphi(x, y) \eta(x, y) \, d\gamma_1^n(x, y) \\
&= \iint_G \int_{\mathbb{R}^d} (K(y, z) - K(x, z)) \varphi(x, y) \eta(x, y) \, d(\rho^n \otimes \gamma_1^n)(z, x, y) \\
&= \iint_{\text{supp } \varphi} \int_{\mathbb{R}^d \cap B_R} (K(y, z) - K(x, z)) \varphi(x, y) \eta(x, y) \, d(\rho^n \otimes \gamma_1^n)(z, x, y) \\
&\quad + \iint_{\text{supp } \varphi} \int_{\mathbb{R}^d \setminus B_R} (K(y, z) - K(x, z)) \varphi(x, y) \eta(x, y) \, d(\rho^n \otimes \gamma_1^n)(z, x, y).
\end{align*}
$$

The last integral is actually vanishing as $R \to \infty$ since (K3), (A1) and Prokhorov's
Theorem give

$$
\begin{align*}
& \left| \iint_{\text{supp } \varphi} \int_{\mathbb{R}^d \setminus B_R} (K(y,z) - K(x,z)) \varphi(x,y) \eta(x,y) \, d(\rho^n \otimes \gamma_1^n)(z,x,y) \right| \\
& \le \frac{L \|\varphi\|_{\infty} \rho^n(\mathbb{R}^d \setminus B_R)}{\inf_{\text{supp } \varphi} (|x-y| \vee |x-y|^2)} \\
& \qquad \times \iint_{\text{supp } \varphi} (|x-y|^2 \vee |x-y|^4) \eta(x,y) \, d\mu^n(y) \, d\rho^n(x) \\
& \le \frac{L C_{\eta} \|\varphi\|_{\infty} \rho^n(\mathbb{R}^d \setminus B_R)}{\inf_{\text{supp } \varphi} (|x-y| \vee |x-y|^2)} \xrightarrow{R \to \infty} 0.
\end{align*}
$$

The function $(z, x, y) \mapsto (K(y, z) - K(x, z))\varphi(x, y)\eta(x, y)$ is continuous and bounded on $(\mathbb{R}^d \cap B_R) \times G$ thanks to Assumption (W). In addition, we note that $(\rho^n \otimes \gamma_1^n)_n$ narrowly converges to $\rho \otimes \gamma_1$ in $\mathcal{P}(\mathbb{R}^d) \times \mathcal{M}^+(G)$. Therefore, we obtain for any $R > 0$ the convergence

$$
\begin{align*}
& \lim_{n \to \infty} \iint_{\text{supp } \varphi} \int_{\mathbb{R}^d \cap B_R} (K(y,z) - K(x,z)) \varphi(x,y) \eta(x,y) \, d(\rho^n \otimes \gamma_1^n)(z,x,y) \\
& = \iint_{\text{supp } \varphi} \int_{\mathbb{R}^d \cap B_R} (K(y,z) - K(x,z)) \varphi(x,y) \eta(x,y) \, d(\rho \otimes \gamma_1)(z,x,y).
\end{align*}
$$

By sending $R \to \infty$, we obtain

$$
\begin{align*}
& \lim_{n \to \infty} \iint_G u^n(x, y) \varphi(x, y) \eta(x, y) \, d\gamma_1^n(x, y) \\
&= \iint_G \iint_{\mathbb{R}^d} (K(y, z) - K(x, z)) \varphi(x, y) \eta(x, y) \, d(\rho \otimes \gamma_1)(z, x, y) \\
&= \iint_G u(x, y) \varphi(x, y) \eta(x, y) \, d\gamma_1(x, y).
\end{align*}
$$
Thus, $u^n$ converges weakly to $u$ as $n \to \infty$ in the sense of [2, Definition 5.4.3]. By [2, Theorem 5.4.4 (ii)] we therefore conclude that

$$ \liminf_{n \to \infty} \mathcal{D}(\mu^n; \rho^n) = \liminf_{n \to \infty} \iint_G g(u^n(x, y))\eta(x, y) \, d\rho^n(x) \, d\mu^n(y) \\ \geq \iint_G g(u(x, y))\eta(x, y) \, d\rho(x) \, d\mu(y) = \mathcal{D}(\mu; \rho), $$

which is the desired result.

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

Let us also prove the compactness and narrow lower semicontinuity of the De Giorgi functional.

**Lemma 3.13.** (Compactness and lower semicontinuity of the De Giorgi functional) Let $(\mu^n)_n \subset \mathcal{M}^+(\mathbb{R}^d)$ and suppose that $(\mu^n)_n$ narrowly converges to $\mu$. Assume that the base measures $\mu^n$ and $\mu$ satisfy (A1) and (A2) uniformly in $n$, and let $K$ satisfy (K1)-(K3). Let moreover $(\rho^n)_n$ be a sequence so that $\rho^n \in AC([0, T]; (\mathcal{P}_2(\mathbb{R}^d), T_{\mu^n}))$ for all $n \in \mathbb{N}$ with $\sup_{n \in \mathbb{N}} M_2(\rho_0^n) < \infty$ and $\sup_{n \in \mathbb{N}} \mathcal{G}_T(\mu^n; \rho^n) < \infty$. Then, up to a subsequence, $\rho_t^n \to \rho_t$ as $n \to \infty$ for all $t \in [0, T]$ for some $\rho \in AC([0, T]; (\mathcal{P}_2(\mathbb{R}^d), T_{\mu^n}))$ and

$$ \liminf_{n \to \infty} \mathcal{G}_T(\mu^n; \rho^n) \geq \mathcal{G}_T(\mu; \rho). $$

**Proof.** For any $n \in \mathbb{N}$, recall the definition

$$ \mathcal{G}_T(\mu^n; \rho^n) = \mathcal{E}(\rho_T^n) - \mathcal{E}(\rho_0^n) + \frac{1}{2} \int_0^T \mathcal{D}(\mu^n; \rho_t^n) \, dt + \frac{1}{2} \int_0^T |(\rho_t^n)'|_{T_{\mu^n}}^2 \, dt, $$

where we are careful to take the metric derivative of $\rho^n$ with respect to $T_{\mu^n}$ (as given in Definition 2.18). Since the domain of the energy $\mathcal{E}$ is all of $\mathcal{P}_2(\mathbb{R}^d)$ and the local slope $\mathcal{D}$ is non-negative, the bound $\sup_{n \in \mathbb{N}} \mathcal{G}_T(\mu^n; \rho^n) < \infty$ ensures that

$$ \sup_{n \in \mathbb{N}} \int_0^T |(\rho_t^n)'|_{T_{\mu^n}}^2 \, dt < \infty. $$

For all $n \in \mathbb{N}$, since $\rho^n \in AC([0, T]; (\mathcal{P}_2(\mathbb{R}^d), T_{\mu^n}))$, Proposition 2.25 yields the existence of a flux $j^n$ such that $(\rho^n, j^n) \in \text{CE}_T$ and $|(\rho_t^n)'|^2 = \mathcal{A}(\mu^n; \rho_t^n, j_t^n)$ for almost all $t \in [0, T]$. We then get

$$ \sup_{n \in \mathbb{N}} \int_0^T \mathcal{A}(\mu^n; \rho_t^n, j_t^n) \, dt = \sup_{n \in \mathbb{N}} \int_0^T |(\rho_t^n)'|^2_{T_{\mu^n}} \, dt < \infty. $$

By Proposition 2.17, there now exists $(\rho, j) \in \text{CE}_T$ such that, up to subsequences, $\rho_t^n \to \rho_t$ for all $t \in [0, T]$ and $j^n \to j$ as $n \to \infty$, and

$$ \infty > \liminf_{n \to \infty} \int_0^T \mathcal{A}(\mu^n; \rho_t^n, j_t^n) \, dt \geq \int_0^T \mathcal{A}(\mu; \rho_t, j_t) \, dt. $$
By Proposition 2.25, we therefore have $\rho \in AC([0, T]; (P_2(\mathbb{R}^d), T_\mu))$ and $|(\rho_t)'|_{T_\mu}^2 \le \mathcal{A}(\mu; \rho_t, j_t)$ for almost all $t \in [0, T]$, which finally gives

$$ \liminf_{n \to \infty} \int_0^T |(\rho_t^n)'|_{T_\mu^n}^2 dt \ge \int_0^T |\rho_t'|_{T_\mu} dt. \quad (3.27) $$

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

By the narrow continuity of the energy proved in Proposition 3.3, we get

$$ \lim_{n \to \infty} \mathcal{E}(\rho_T^n) = \mathcal{E}(\rho_T) \quad \text{and} \quad \lim_{n \to \infty} \mathcal{E}(\rho_0^n) = \mathcal{E}(\rho_0). \quad (3.28) $$

Furthermore, by Fatou's lemma and the narrow lower semicontinuity of the local slope shown in Lemma 3.12, we have

$$ \liminf_{n \to \infty} \int_0^T \mathcal{D}(\mu^n; \rho_t^n) dt \ge \int_0^T \mathcal{D}(\mu; \rho_t) dt. \quad (3.29) $$

Gathering (3.27), (3.28) and (3.29), we finally obtain

$$ \begin{aligned} \liminf_{n \to \infty} \mathcal{G}_T(\mu^n; \rho^n) &\ge \mathcal{E}(\rho_T) - \mathcal{E}(\rho_0) + \frac{1}{2} \int_0^T \mathcal{D}(\mu; \rho_t) dt + \frac{1}{2} \int_0^T |\rho_t'|_{T_\mu}^2 dt \\ &= \mathcal{G}_T(\mu; \rho), \end{aligned} $$

which ends the proof.

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

We now get our stability result.

**Theorem 3.14.** (Stability of gradient flows) Let $(\mu^n)_n \subset \mathcal{M}^+(\mathbb{R}^d)$ and suppose that $(\mu^n)_n$ narrowly converges to $\mu$. Assume that the base measures $\mu^n$ and $\mu$ satisfy (A1) and (A2) uniformly in $n$, and let the interaction potential $K$ satisfy (K1)–(K3). Suppose that $\rho^n$ is a gradient flow of $\mathcal{E}$ with respect to $\mu^n$ for all $n \in \mathbb{N}$, that is,

$$ \mathcal{G}_T(\mu^n; \rho^n) = 0 \quad \text{for all } n \in \mathbb{N}, $$

such that $(\rho_0^n)_n$ satisfies $\sup_{n \in \mathbb{N}} M_2(\rho_0^n) < \infty$ and $\rho_t^n \to \rho_t$ as $n \to \infty$ for all $t \in [0, T]$ for some curve $(\rho_t)_{t \in [0, T]} \subset \mathcal{P}_2(\mathbb{R}^d)$. Then, $\rho \in AC([0, T]; (P_2(\mathbb{R}^d), T_\mu))$ and $\rho$ is a gradient flow of $\mathcal{E}$ with respect to $\mu$, that is,

$$ \mathcal{G}_T(\mu; \rho) = 0. $$

**Proof.** By Lemma 3.13 we directly obtain that $\rho \in AC([0, T]; (P_2(\mathbb{R}^d), T_\mu))$ and, up to a subsequence,

$$ 0 = \liminf_{n \to \infty} \mathcal{G}_T(\mu^n; \rho^n) \ge \mathcal{G}(\mu; \rho). $$

Finally, since $\mathcal{G}_T(\mu; \rho) \ge 0$ by Young's inequality and Corollary 3.11, we obtain $\mathcal{G}_T(\mu; \rho) = 0$. □

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

Note that, via Theorem 3.9, the above theorem also shows stability of weak solutions to (NL$^2$IE). Typically, in Theorem 3.14, $(\mu^n)_n$ is a sequence of atomic measures used to approximate, or sample, the support of $\mu$. Indeed, we now use this approach to show the existence of weak solutions to the nonlocal nonlocal-interaction equation.
**Theorem 3.15.** (Existence of weak solutions) Let $K$ be an interaction potential satisfying Assumptions (K1)–(K3). Suppose that $\mu \in \mathcal{M}^+(\mathbb{R}^d)$ is finite, i.e., $\mu(\mathbb{R}^d) < \infty$, and satisfies (A2). Assume furthermore that for some $C'_\eta > 0$ it holds that

$$ \sup_{(x,y) \in G \cap \text{supp } \mu \otimes \mu} \left( |x - y|^2 \lor |x - y|^4 \right) \eta(x, y) \leq C'_\eta. \quad (3.30) $$

Consider $\rho_0 \in \mathcal{P}_2(\mathbb{R}^d)$ which is $\mu$-absolutely continuous. Then there exists a weakly continuous curve $\rho: [0, T] \to \mathcal{P}(\mathbb{R}^d)$ such that $\text{supp } \rho_t \subseteq \text{supp } \mu$ for all $t \in [0, T]$, which is a weak solution of (NL$^2$IE) and satisfies the initial condition $\rho(0) = \rho_0$.

**Proof.** Let $(\mu^n)_n \subset \mathcal{M}^+(\mathbb{R}^d)$ be a sequence of atomic measures such that $(\mu^n)_n$ converges narrowly to $\mu$. Moreover, assume that $\mu^n$ has finitely many atoms and $\mu^n(\mathbb{R}^d) \leq \mu(\mathbb{R}^d)$ and $\text{supp } \mu^n \subseteq \text{supp } \mu$ for all $n \in \mathbb{N}$. Let $\hat{\mu}^n$ be the normalization of $\mu^n$ which has the same total mass as $\mu$, that is,

$$ \hat{\mu}^n = \frac{\mu(\mathbb{R}^d)}{\mu^n(\mathbb{R}^d)} \mu^n, $$

and let $\pi^n$ be optimal transportation plan between $\mu$ and $\hat{\mu}^n$ for the quadratic cost. Let $\rho_0^n$ be the second marginal of $\tilde{\rho}_0\pi^n$, where $\tilde{\rho}_0$ is the density of the measure $\rho_0$ with respect to $\mu$; namely, let $\rho_0^n(A) = \int_{\mathbb{R}^d \times A} \tilde{\rho}_0(x) \, d\pi^n(x, y)$ for any Borel set $A \subset \mathbb{R}^d$. Note that $\rho_0^n(\mathbb{R}^d) = \rho_0(\mathbb{R}^d)$ and $\rho_0^n \ll \mu^n$ for all $n \in \mathbb{N}$, and that, since $\tilde{\rho}_0\pi^n$ is a transport plan between $\rho_0$ and $\rho_0^n$, $\rho_0^n \to \rho_0$ as $n \to \infty$.

Thanks to Assumption (3.30), it holds, for all $n \in \mathbb{N}$, that

$$ \mu - \text{ess sup}_{x \in \mathbb{R}^d} \int (|x - y|^2 \lor |x - y|^4) \eta(x, y) \, d\mu^n(y) \leq \mu^n(\mathbb{R}^d) C'_\eta \\ \leq \mu(\mathbb{R}^d) C'_\eta. \quad (3.31) $$

Since, by construction $\rho_0^n \ll \mu^n$, we have $\text{supp } \rho_0^n \subseteq \text{supp } \mu^n \subseteq \text{supp } \mu$. This nested support property is, thanks to Proposition 2.28, preserved in time, so that $\text{supp } \rho_t^n \subseteq \text{supp } \mu$ for all $t \in [0, T]$ and $n \in \mathbb{N}$. For this reason, (3.31) can be used, under the stated support restriction on $\rho_0$, instead of Assumption (A1) uniformly in $n$ when calling Lemma 3.13 and Theorem 3.14 later in this proof. Since $\mu^n$ consists of finitely many atoms and $\mu$ satisfies (A2), the family $(\mu_n)_n$ satisfies (A2) uniformly in $n$.

By Remark 1.1, we know that the ODE system (1.2)–(1.4) admits a unique solution for all $n \in \mathbb{N}$. It can be easily checked that this solution, which we denote by $\rho^n$, is a weak solution to (NL$^2$IE) with respect to $\mu^n$ starting from $\rho_0^n$, according to Definition 3.1. By Theorem 3.9, we then get that $\rho^n$ is a gradient flow of $\mathcal{E}$ with respect to $\mu$ starting from $\rho_0^n$ for all $n \in \mathbb{N}$.

Combining the compactness part of Lemma 3.13 and the stability from Theorem 3.14, we get that, up to a subsequence, $\rho_t^n \to \rho_t$ as $n \to \infty$ for all $t \in [0, T]$, where $\rho \in \text{AC}([0, T]; (\mathcal{P}_2(\mathbb{R}^d), \mathcal{T}_\mu))$ is a gradient flow of $\mathcal{E}$ with respect to $\mu$ starting from $\rho_0$. Theorem 3.9 finally shows that $\rho$ is a weak solution to (NL$^2$IE) with respect to $\mu$ starting from $\rho_0$. □
**Remark 3.16.** Assumption (3.30) is only needed to arrive at an atomic approximation sequence $(\mu^n)_n$ of $\mu$ such that Assumptions (A1) and (A2) hold uniformly in $n$. On a case-by-case basis, one could drop (3.30) and try to construct the sequence $(\mu^n)_n$ explicitly in such a way as to satisfy both assumptions uniformly in $n$.

**Remark 3.17.** We conclude the section by remarking on the relevance of the Theorem 3.14 to the setting of machine learning. Namely, there $\mu$ is the measure modeling the true data distribution, which can be assumed to be compact. Let $(x_i)_i$ be a sequence of i.i.d. samples of $\mu$ and let $\mu^n = \frac{1}{n} \sum_{i=1}^n \delta_{x_i}$ be the empirical measure of the first $n$ sample points. Assume $(\rho^n)_n$ is a narrowly converging sequence of probability measures such that $\supp \rho^n \subseteq \{x_1, \dots, x_n\}$ for all $n \in \mathbb{N}$, and denote by $\rho$ its limit. Assume that $\eta$ is an edge weight kernel such that $\mu$ and $\eta$ satisfy (A2) and (3.30). Let $K$ be an interaction kernel satisfying (K2) and (K3). Finally, let $(\tilde{\rho}^n)_n$ be the sequence of solutions of (NL$^2$IE) in the sense of Definition 3.1 such that $\tilde{\rho}_0^n = \rho^n$ for all $n \in \mathbb{N}$. Then, by Lemma 3.13, the sequence $(\tilde{\rho}_t^n)_n$ narrowly converges along a subsequence for all $t \in [0, T]$, and furthermore, by Theorem 3.15, any curve $(\tilde{\rho}_t)_{t \in [0,T]}$ of subsequential limits yields a solution $\tilde{\rho}$ of (NL$^2$IE) with initial condition $\rho$.

### 3.5. Discussion of the Local Limit

Here we discuss at a formal level the connection between the nonlocal nonlocal-interaction equation and its limit as the graph structure localizes. We first present a very formal justification as to why we expect the solutions of (NL$^2$IE) to converge to the solutions of a nonlocal-interaction equation as the localizing parameter $\varepsilon \to 0^+$, i.e., as the edge-weight function $\eta = \eta_\varepsilon$ localizes. We conclude this section with an example that cautions that the formal argument cannot be justified in full generality. Proving the convergence of (NL$^2$IE) in the limit $\varepsilon \to 0^+$, under appropriate conditions, remains an intriguing open problem.

Take $\mu = \text{Leb}(\mathbb{R}^d)$ and choose $\eta_\varepsilon$ given by (2.2). Consider a smooth interaction potential $K: \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}$ and a compactly supported initial condition $\rho_0$ which has a continuous density with respect to $\mu$. Let $\rho^\varepsilon$ be the solution of (NL$^2$IE) starting from $\rho_0$ for the edge weight function $\eta_\varepsilon$. Assume that $\rho_t^\varepsilon$ is absolutely continuous with respect to $\mu$ for all $t$. In the following we drop the $t$-dependence of $\rho^\varepsilon$ for brevity. From (NL$^2$IE), by adding and subtracting $\rho^\varepsilon(x) \int_{\mathbb{R}^d} (\bar{\nabla} K * \rho^\varepsilon(x, y))_+ \eta_\varepsilon(x, y) \, dy$, it follows that

$$ \begin{aligned} \partial_t \rho^\varepsilon(x) = & -\rho^\varepsilon(x) \int_{\mathbb{R}^d} \bar{\nabla} K * \rho^\varepsilon(x, y) \eta_\varepsilon(x, y) \, dy \\ & - \int_{\mathbb{R}^d} \bar{\nabla} \rho^\varepsilon(x, y) (\bar{\nabla} K * \rho^\varepsilon(x, y))_+ \eta_\varepsilon(x, y) \, dy. \end{aligned} $$

Then, for almost all $x \in \mathbb{R}^d$ we have

$$ \int_{\mathbb{R}^d} \bar{\nabla} K * \rho^\varepsilon(x, y) \eta_\varepsilon(x, y) \, dy $$
$$ 
\begin{aligned}
&= \frac{2(2+d)}{\varepsilon^2} \int_{\mathbb{R}^d} (K * \rho^\varepsilon(y) - K * \rho^\varepsilon(x)) \frac{\chi_{B_\varepsilon}(x)(y)}{|B_\varepsilon|} \mathrm{d}y \\
&= \frac{2(2+d)}{\varepsilon^2} \left( \frac{1}{|B_\varepsilon|} \int_{B_\varepsilon(x)} K * \rho^\varepsilon(y) \mathrm{d}y - K * \rho^\varepsilon(x) \right).
\end{aligned}
$$

A standard calculation, using a second-order Taylor expansion, shows that the right-hand side approximates $\Delta K * \rho^\varepsilon(x)$ when $\varepsilon$ is small, provided that derivatives of $\rho^\varepsilon$ remain uniformly bounded.

Similarly, by Taylor expanding $\overline{\nabla} \rho^\varepsilon$ and $\overline{\nabla} K * \rho^\varepsilon$ to first order and changing variable over the unit sphere while carefully tracking the positive part, one gets

$$ \int_{\mathbb{R}^d} \overline{\nabla} \rho^\varepsilon(x, y) (\overline{\nabla} K * \rho^\varepsilon(x, y))_+ \eta_\varepsilon(x, y) \mathrm{d}y \approx \overline{\nabla} \rho^\varepsilon(x) \cdot \overline{\nabla} K * \rho^\varepsilon(x) $$
for small $\varepsilon$.

Combining the expressions above yields

$$ 
\begin{aligned}
\partial_t \rho^\varepsilon(x) &\approx -\rho^\varepsilon(x) \Delta K * \rho^\varepsilon(x) - \overline{\nabla} \rho^\varepsilon(x) \cdot \overline{\nabla} K * \rho^\varepsilon(x) \\
&= -\overline{\nabla} \cdot (\rho^\varepsilon \overline{\nabla} K * \rho^\varepsilon)(x).
\end{aligned}
$$

This suggests that if $\rho^\varepsilon$ converge as $\varepsilon \to 0^+$, then the limiting $\rho$ is a solution of the standard nonlocal interaction equation (3.1). A possible way to attack the local limit within the variational framework is via a stability statement similar to that of Theorem 3.14, but now with respect to the family $(\eta_\varepsilon)_{\varepsilon>0}$ in the limit $\varepsilon \to 0^+$. The next remark indicates that this will require further regularity assumptions on the interaction kernel $K$.

**Remark 3.18.** We present an example that indicates that, in certain situations, solutions of (NL$^2$IE) cannot be expected to converge to solutions of (3.1) as the interaction kernel $\eta_\varepsilon$ becomes more concentrated. Namely, consider $d = 1$, $\Omega = (-2, 2)$ and $\mu = \text{Leb}(\Omega)$. Let $K(x, y) = 1 - e^{-|x-y|}$ for all $x, y \in \Omega$ and $\eta$ be a smooth, even function, positive on $(-0.2, 0.2)$ and zero otherwise. Consider $\rho_0 = \frac{1}{2}(\delta_{-1} + \delta_1)$. It is straightforward to verify that $\rho_t = \rho_0$ for all $t \in [0, T]$ yields a weak solution of (NL$^2$IE) for all $\varepsilon > 0$. In particular, note that the corresponding velocity field satisfies $v(-1, y) = -(K * \rho_0(y) - K * \rho_0(-1)) \le 0$ for all $y \in (-1.2, -0.8)$, and thus the flux from $x = -1$ remains zero, and analogously from $x = 1$. Therefore, one cannot expect the weak solutions for the interaction potential $K$ to converge to weak solutions of (3.1) as $\varepsilon \to 0^+$. We believe that, for these particular kernel $K$ and edge weights $\eta$, the problem persists for strong solutions for initial data close to $\rho_0$, only that explicit solutions are not available.

*Acknowledgements.* The authors are grateful to José Antonio CARRILLO for several enlightening discussions, to MARK PELETIER for insightful remarks on an earlier version of the paper, and to TRIPHON GEORGIOU for valuable information. The authors would like to thank the American Institute of Mathematics (AIM) for its support through the workshop *Nonlocal differential equations in collective behavior*. AE and AS gratefully acknowledge the support of the Hausdorff Research Institute for Mathematics (Bonn), through the Junior Trimester Program on *Kinetic Theory*. AE acknowledges support by the EU-funded Erasmus Mundus programme "MathMod - Mathematical models in engineering: theory, methods, and applications" at the University of L'Aquila, and by the German Science Foundation (DFG) through
CRC TR 154 "Mathematical Modelling, Simulation and Optimization Using the Example of Gas Networks". AS is supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy EXC 2044 - 390685587, Mathematics Münster: Dynamics–Geometry–Structure, and EXC 2047 - 390685813, Hausdorff Center for Mathematics, as well as the Collaborative Research Center 1060 - 211504053, The Mathematics of Emergent Effects at the Universität Bonn. DS is grateful to NSF for support via Grants DMS 1516677 and DMS 1814991 and via Ki-Net (SF Research Network Grant RNMS 1107444). DS and FSP are grateful to the Center for Nonlinear Analysis of CMU for its support.

**Funding** Open Access funding enabled and organized by Projekt DEAL.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Appendix A: Minkowski Norm of the Underlying Finsler Structure

In this appendix we show that, given $\rho \in \mathcal{P}_2(\mathbb{R}^d)$ and $j \in T_\rho \mathcal{P}_2(\mathbb{R}^d)$, the inner product $g_{\rho,j}$ from Section 3.1 derives from a so-called Minkowski norm, as it should be in the theory of Finsler geometry; see [4, 18, 49].

Let us fix $\rho \in \mathcal{P}_2(\mathbb{R}^d)$ which is absolutely continuous with respect to $\mu$, in accordance with Section 3.1. For $j \in T_\rho \mathcal{P}_2(\mathbb{R}^d)$, we denote $j$ its density with respect to $\mu \otimes \mu$. We show that the function $F_\rho: T_\rho \mathcal{P}_2(\mathbb{R}^d) \to \mathbb{R}$ given by

$$F_{\rho}(j) = \sqrt{\frac{1}{2} \iint_{G} j(x, y)^2 \eta(x, y) \left( \frac{\chi_{\{j>0\}}(x, y)}{\rho(x)} + \frac{\chi_{\{j<0\}}(x, y)}{\rho(y)} \right) d\mu(x) d\mu(y)} \\ \text{for all } j \in T_{\rho} \mathcal{P}_2(\mathbb{R}^d),$$

is a Minkowski norm, that is, it is smooth away from 0, positively 1-homogeneous and, whenever $j$ in nonzero $\eta \mu \otimes \mu$-a.e., its second variation is a symmetric positive definite bilinear form. In fact, we now prove that, for all $j, j_1, j_2 \in T_\rho \mathcal{P}_2(\mathbb{R}^d)$ such that $j$ is nonzero $\eta \mu \otimes \mu$-a.e.:

$$\left. \frac{1}{2} \frac{\partial^2}{\partial t \partial s} \right|_{s=t=0} F_{\rho}^2(j + sj_1 + tj_2) = g_{\rho,j}(j_1, j_2).$$

Indeed, let $j \in T_\rho \mathcal{P}_2(\mathbb{R}^d)$, $s \in \mathbb{R}$ and $j_1 \in T_\rho \mathcal{P}_2(\mathbb{R}^d)$ such that $j + sj_1 \in T_\rho \mathcal{P}_2(\mathbb{R}^d)$. Then,

$$F_{\rho}^{2}(j + sj_{1}) - F_{\rho}^{2}(j) \\ = \frac{1}{2} \iint_{G} (j(x, y) + sj_{1}(x, y))^{2} \eta(x, y) \\ \times \left( \frac{\chi_{\{j + sj_{1} > 0\}}(x, y)}{\rho(x)} + \frac{\chi_{\{j + sj_{1} < 0\}}(x, y)}{\rho(y)} \right) d\mu(x) d\mu(y)$$
$$
\begin{align*}
& -\frac{1}{2} \iint_G j(x, y)^2 \eta(x, y) \left( \frac{\chi_{\{j>0\}}(x, y)}{\rho(x)} + \frac{\chi_{\{j<0\}}(x, y)}{\rho(y)} \right) d\mu(x) d\mu(y) \\
&= \frac{1}{2} \iint_G j(x, y)^2 \eta(x, y) \left( \frac{\chi_{\{j+sj_1>0\}}(x, y)}{\rho(x)} - \frac{\chi_{\{j>0\}}(x, y)}{\rho(x)} \right) d\mu(x) d\mu(y) \\
&\quad + \frac{1}{2} \iint_G j(x, y)^2 \eta(x, y) \left( \frac{\chi_{\{j+sj_1<0\}}(x, y)}{\rho(y)} - \frac{\chi_{\{j<0\}}(x, y)}{\rho(y)} \right) d\mu(x) d\mu(y) \\
&\quad + \frac{s^2}{2} \iint_G j_1(x, y)^2 \eta(x, y) \\
&\quad \quad \times \left( \frac{\chi_{\{j+sj_1>0\}}(x, y)}{\rho(x)} + \frac{\chi_{\{j+sj_1<0\}}(x, y)}{\rho(y)} \right) d\mu(x) d\mu(y) \\
&\quad + s \iint_G j(x, y) j_1(x, y) \eta(x, y) \\
&\quad \quad \times \left( \frac{\chi_{\{j+sj_1>0\}}(x, y)}{\rho(x)} + \frac{\chi_{\{j+sj_1<0\}}(x, y)}{\rho(y)} \right) d\mu(x) d\mu(y) \\
&=: I_1 + I_2 + I_3 + I_4.
\end{align*}
$$

Note that

$$
\begin{align*}
I_1 &= \frac{1}{2} \iint_{\{0 < j \le -sj_1\}} j(x, y)^2 \eta(x, y) \\
    &\quad \times \left( \frac{\chi_{\{j+sj_1>0\}}(x, y)}{\rho(x)} - \frac{\chi_{\{j>0\}}(x, y)}{\rho(x)} \right) d\mu(x) d\mu(y) \\
    &\quad + \frac{1}{2} \iint_{\{-sj_1 < j \le 0\}} j(x, y)^2 \eta(x, y) \\
    &\quad \times \left( \frac{\chi_{\{j+sj_1>0\}}(x, y)}{\rho(x)} - \frac{\chi_{\{j>0\}}(x, y)}{\rho(x)} \right) d\mu(x) d\mu(y) \\
    &= -\frac{1}{2} \iint_{\{0 < j \le -sj_1\}} \frac{j(x, y)^2 \eta(x, y)}{\rho(x)} d\mu(x) d\mu(y) \\
    &\quad + \frac{1}{2} \iint_{\{-sj_1 < j < 0\}} \frac{j(x, y)^2 \eta(x, y)}{\rho(x)} d\mu(x) d\mu(y).
\end{align*}
$$

Therefore,

$$
|I_1| \leq s^2 \iint_G \frac{j_1(x, y)^2 \eta(x, y)}{\rho(x)} d\mu(x) d\mu(y).
$$

Similarly, one gets

$$
|I_2| \leq s^2 \iint_G \frac{j_1(x, y)^2 \eta(x, y)}{\rho(y)} d\mu(x) d\mu(y).
$$

We also have that

$$
|I_3| \le \frac{s^2}{2} \iint_G \frac{j_1(x, y)\eta(x, y)}{\rho(x)} d\mu(x) d\mu(y) \\
+ \frac{s^2}{2} \iint_G \frac{j_1(x, y)\eta(x, y)}{\rho(y)} d\mu(x) d\mu(y).
$$

Using Lebesgue's dominated convergence theorem, one gets, moreover, that

$$
\lim_{s \to 0} \frac{I_4}{s} = \iint_G j_1(x, y) j(x, y) \eta(x, y)
$$
$$ \times \left( \frac{\chi_{\{j>0\}}(x, y)}{\rho(x)} + \frac{\chi_{\{j<0\}}(x, y)}{\rho(y)} \right) d\mu(x) d\mu(y). $$

We thus overall get

$$ \lim_{s \to 0} \frac{F_{\rho}^{2}(j + s j_{1}) - F_{\rho}^{2}(j)}{s} = \iint_{G} j_{1}(x, y) j(x, y) \eta(x, y) \times \left( \frac{\chi_{\{j>0\}}(x, y)}{\rho(x)} + \frac{\chi_{\{j<0\}}(x, y)}{\rho(y)} \right) d\mu(x) d\mu(y), $$

which shows that

$$ \left. \frac{1}{2} \frac{d}{ds} \right|_{s=0} F_{\rho}^{2}(j + s j_{1}) = \frac{1}{2} \iint_{G} j_{1}(x, y) j(x, y) \eta(x, y) \times \left( \frac{\chi_{\{j>0\}}(x, y)}{\rho(x)} + \frac{\chi_{\{j<0\}}(x, y)}{\rho(y)} \right) d\mu(x) d\mu(y) \\ = g_{\rho, j}(j_{1}, j). \qquad (A.1) $$

Note that this equality was used in Section 3.1 to determine that grad⁻ $\mathcal{E}(\rho)$ is indeed the direction of steepest descent from $\rho$, i.e., to get (3.9). Computing now a further derivative in direction $j_2 \in T_\rho P_2(\mathbb{R}^d)$ and using similar boundedness arguments, we get

$$ \left. \frac{1}{2} \frac{\partial^2}{\partial t \partial s} \right|_{s=t=0} F_{\rho}^{2}(j + s j_{1} + t j_{2}) \\ = \frac{1}{2} \iint_{G} j_{1}(x, y) j_{2}(x, y) \eta(x, y) \times \left( \frac{\chi_{\{j>0\}}(x, y)}{\rho(x)} + \frac{\chi_{\{j<0\}}(x, y)}{\rho(y)} \right) d\mu(x) d\mu(y) \\ + \frac{1}{2} \iint_{\{j=0\}} j_{1}(x, y) j_{2}(x, y) \eta(x, y) \times \left( \frac{\chi_{\{j>0\}}(x, y)}{\rho(x)} + \frac{\chi_{\{j<0\}}(x, y)}{\rho(y)} \right) d\mu(x) d\mu(y) \\ = g_{\rho, j}(j_{1}, j_{2}) + \frac{1}{2} \iint_{\{j=0\}} j_{1}(x, y) j_{2}(x, y) \eta(x, y) \times \left( \frac{\chi_{\{j>0\}}(x, y)}{\rho(x)} + \frac{\chi_{\{j<0\}}(x, y)}{\rho(y)} \right) d\mu(x) d\mu(y). $$

Note that the presence of the integral over the set $\{j = 0\}$ comes from the fact that $j$ is not a multiplicative function of the integrand anymore (as it was the case for the first derivative), so that the set of points where $j = 0$ has to be considered. Assuming then that $\eta \mu \otimes \mu$-a.e. we have $j \neq 0$ we obtain that this integral over $\{j = 0\}$ is equal to zero, which yields the claim.

## References

1. AGUEH, M.: Finsler structure in the $p$-Wasserstein space and gradient flows. C. R. Math. Acad. Sci. Paris 350(1–2), 35–40, 2012

2. AMBROSIO, L., GIGLI, N., SAVARÉ, G.: Gradient Flows in Metric Spaces and in the Space of Probability Measures. Lectures in Mathematics ETH Zurich. Birkhäuser Verlag, Basel 2008
3. BALAGUÉ, D., CARRILLO, J.A., LAURENT, T., RAOUL, G.: Dimensionality of local minimizers of the interaction energy. *Arch. Ration. Mech. Anal.* **209**(3), 1055–1088, 2013

4. BAO, D., CHERN, S.-S., SHEN, Z.: *An Introduction to Riemann-Finsler Geometry*. Springer, New York 2000

5. BELKIN, M., NIYOGI, P.: Laplacian eigenmaps for dimensionality reduction and data representation. *Neural Comput.* **15**, 1373–1396, 2002

6. BENAMOU, J.-D., BRENIER, Y.: A computational fluid mechanics solution to the Monge–Kantorovich mass transfer problem. *Numer. Math.* **84**(3), 375–393, 2000

7. BERTOZZI, A.L., CARRILLO, J.A., LAURENT, T.: Blow-up in multidimensional aggregation equations with mildly singular interaction kernels. *Nonlinearity* **22**(3), 683–710, 2009

8. BUTTAZZO, G.: *Semicontinuity, relaxation and integral representation in the calculus of variations, volume 207 of Pitman Research Notes in Mathematics Series*. Longman Scientific & Technical, Harlow; copublished in the USA with John Wiley & Sons Inc, New York, 1989

9. CANCÈS, C., GALLOUËT, T.O., TODESCHI, G.: A variational finite volume scheme for Wasserstein gradient flows. Preprint arXiv:1907.08305, 2019

10. CARRILLO, J.A., DI FRANCESCO, M., FIGALLI, A., LAURENT, T., SLEPCEV, D.: Global-in-time weak measure solutions and finite-time aggregation for nonlocal interaction equations. *Duke Math. J.* **156**(2), 229–271, 2011

11. CARRILLO, J.A., McCANN, R.J., VILLANI, C.: Kinetic equilibration rates for granular media and related equations: entropy dissipation and mass transportation estimates. *Rev. Mat. Iberoam.* **19**(3), 971–1018, 2003

12. CHAUDHURI, K., DASGUPTA, S., KPOTUFE, S., VON LUXBURG, U.: Consistent procedures for cluster tree estimation and pruning. *IEEE Trans. Inform. Theory* **60**(12), 7900–7912, 2014

13. CHEN, Y., GEORGIOU, T.T., TANNENBAUM, A.: Vector-valued optimal mass transport. *SIAM J. Appl. Math.* **78**(3), 1682–1696, 2018

14. CHOW, S.-N., HUANG, W., LI, Y., ZHOU, H.: Fokker-Planck equations for a free energy functional or Markov process on a graph. *Arch. Ration. Mech. Anal.* **203**(3), 969–1008, 2012

15. CHOW, S.-N., LI, W., ZHOU, H.: Entropy dissipation of Fokker-Planck equations on graphs. *Discrete Contin. Dyn. Syst.* **38**(10), 4929–4950, 2018

16. COIFMAN, R.R., LAFON, S.: Diffusion maps. *Appl. Comput. Harmon. Anal.* **21**(1), 5–30, 2006

17. COMANICIU, D., MEER, P.: Mean shift: a robust approach toward feature space analysis. *IEEE Trans. Pattern Anal. Mach. Intell.* **24**(5), 603–619, 2002

18. DAHL, M.: A brief introduction to Finsler geometry. 2006. This work is based on my licentiate thesis, Propagation of Gaussian beams using Riemann-Finsler geometry, Helsinki University of technology, 2006, instructed by doctor Kirsi Peltonen and supervised by professor Erkki Somersalo

19. DELARUE, F., LAGOUTIÈRE, F.: Probabilistic analysis of the upwind scheme for transport equations. *Arch. Ration. Mech. Anal.* **199**(1), 229–268, 2011

20. DELARUE, F., LAGOUTIÈRE, F., VAUCHELET, N.: Convergence order of upwind type schemes for transport equations with discontinuous coefficients. *J. Math. Pures Appl.*, 2017 (to appear)

21. DISSER, K., LIERO, M.: On gradient structures for Markov chains and the passage to Wasserstein gradient flows. *Netw. Heterog. Media* **10**(2), 233–253, 2015

22. DOLBEAULT, J., NAZARET, B., SAVARÉ, G.: A new class of transport distances between measures. *Calc. Var. Partial Differ. Equ.* **34**(2), 193–231, 2009

23. ERBAR, M.: Gradient flows of the entropy for jump processes. *Ann. Inst. Henri Poincaré Probab. Stat.* **50**(3), 920–945, 2014

24. ERBAR, M.: A gradient flow approach to the Boltzmann equation. Preprint arXiv:1603.0540, 2018
25. ERBAR, M., FATHI, M., LASCHOS, V., SCHLICHTING, A.: Gradient flow structure for McKean-Vlasov equations on discrete spaces. *Discrete Contin. Dyn. Syst.* **36**(12), 6799–6833, 2016

26. ERBAR, M., FATHI, M., SCHLICHTING, A.: Entropic curvature and convergence to equi-librium for mean-field dynamics on discrete spaces. *ALEA Lat. Am. J. Probab. Math. Stat.* **17**(1), 445–471, 2020

27. ERBAR, M., MAAS, J.: Gradient flow structures for discrete porous medium equations. *Discrete Contin. Dyn. Syst.* **34**(4), 1355–1374, 2014

28. ERBAR, M., MAAS, J., WIRTH, M.: On the geometry of geodesics in discrete optimal transport. *Calc. Var. Partial Differ. Equ.* **58**(1), 19, 2019

29. EYMARD, R., GALLOUËT, T., HERBIN, R.: Finite volume methods. In: Handbook of Numerical Analysis, number Part 3, pp. 713–1018. Elsevier, 2000

30. GANGBO, W., LI, W., MOU, C.: Geodesics of minimal length in the set of probability measures on graphs. *ESAIM Control Optim. Calc. Var.* **25**, 78, 2019

31. GARCÍA TRILLOS, N., SLEPCEV, D.: Continuum limit of total variation on point clouds. *Arch. Ration. Mech. Anal.* **220**(1), 193–241, 2016

32. GARCÍA TRILLOS, N., SLEPCEV, D.: A variational approach to the consistency of spectral clustering. *Appl. Comput. Harmon. Anal.* **45**(2), 239–281, 2018

33. JORDAN, R., KINDERLEHRER, D., OTTO, F.: The variational formulation of the Fokker-Planck equation. *SIAM J. Math. Anal.* **29**(1), 1–17, 1998

34. KANNAN, R., VEMPALA, S., VETTA, A.: On clusterings: good, bad and spectral. *J. ACM* **51**(3), 497–515, 2004

35. KOLOKOLNIKOV, T., SUN, H., UMINSKY, D., BERTOZZI, A.: Stability of ring patterns arising from two-dimensional particle interactions. *Phys. Rev. E Stat. Nonlinear Soft Matter Phys.* **84** 1 Pt 2, 015203, 2011

36. MAAS, J.: Gradient flows of the entropy for finite Markov chains. *J. Funct. Anal.* **261**(8), 2250–2292, 2011

37. MIELKE, A.: A gradient structure for reaction-diffusion systems and for energy-drift-diffusion systems. *Nonlinearity* **24**(4), 1329–1346, 2011

38. MIELKE, A.: Geodesic convexity of the relative entropy in reversible Markov chains. *Calc. Var. Partial Differ. Equ.* **48**(1–2), 1–31, 2013

39. NATALE, A., TODESCHI, G.: TPFA finite volume approximation of Wasserstein gradient flows. In: Finite Volumes for Complex Applications IX - Methods, Theoretical Aspects, Examples, pp. 193–201. Springer 2020

40. NG, A. Y., JORDAN, M. I., WEISS, Y.: On spectral clustering: analysis and an algorithm. In: Advances in Neural Information Processing Systems, pp. 849–856. MIT Press 2001

41. OHTA, S.-I., STURM, K.-T.: Heat flow on Finsler manifolds. *Commun. Pure Appl. Math.* **62**(10), 1386–1433, 2009

42. OHTA, S.-I., STURM, K.-T.: Non-contraction of heat flow on Minkowski spaces. *Arch. Ration. Mech. Anal.* **204**(3), 917–944, 2012

43. PELETIER, M. A., ROSSI, R., SAVARÉ, G., TSE, O.: Jump processes as Generalized Gradient Flows. Preprint arXiv:2006.10624, 2020

44. SCHARFETTER, D. L., GUMMEL, H. K.: Large-signal analysis of a silicon read diode oscillator. *IEEE Trans. Electron Devices* **16**(1), 64–77, 1969

45. SCHLICHTING, A., SEIS, C.: Convergence rates for upwind schemes with rough coefficients. *SIAM J. Numer. Anal.* **55**(2), 812–840, 2017

46. SCHLICHTING, A., SEIS, C.: Analysis of the implicit upwind finite volume scheme with rough coefficients. *Numer. Math.* **139**(1), 155–186, 2018

47. SCHLICHTING, A., SEIS, C.: The Scharfetter-Gummel scheme for aggregation-diffusion equations. Preprint arXiv:2004.13981, 2020

48. SERFATY, S.: Gamma-convergence of gradient flows on Hilbert and metric spaces and applications. *Discrete Contin. Dyn. Syst.* **31**(4), 1427–1451, 2011

49. SHEN, Z.: *Lectures on Finsler Geometry*. World Scientific, Singapore 2001

50. VILLANI, C.: Topics in Optimal Transportation, Volume 58 of Graduate Studies in Mathematics. American Mathematical Society, Providence 2003
A. ESPOSITO
Department Mathematik,
Friedrich-Alexander-Universität Erlangen-Nürnberg,
Cauerstrasse 11,
91058 Erlangen
Germany.
e-mail: antonio.esposito@fau.de
and

F. S. PATACCHINI · D. SLEPCEV
Department of Mathematical Sciences,
Carnegie Mellon University,
Pittsburgh
PA
15213 USA.
e-mail: slepcev@math.cmu.edu
and

F. S. PATACCHINI
IFP Energies nouvelles,
1 et 4 avenue de Bois-Préau,
92852 Rueil-Malmaison
France.
e-mail: fpatacch@math.cmu.edu
and
e-mail: francesco.patacchini@ifpen.fr
and

A. SCHLICHTING
Institute for Analysis and Numerics,
University of Münster,
Orléans-Ring 10,
48149 Münster
Germany.
e-mail: schlichting@iam.uni-bonn.de
and

Institute for Applied Mathematics,
University of Bonn,
Endenicher Allee 60,
53115 Bonn
Germany.
and

*Present address:* Institute for Analysis and Numerics,
University of Münster,
Orléans-Ring 10,
48149 Münster
Germany.
e-mail: a.schlichting@uni-muenster.de